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ABSTRACT 



We investigate the time evolution of dust properties, molecular hydrogen (H2) 
contents, and star formation histories in galaxies by using our original chemodynami- 
cal simulations. The simulations include the formation of dust in the stellar winds of 
supernovae (SNe) and asymptotic giant branch (AGB) stars, the growth and destruc- 
tion processes of dust in the interstellar medium (ISM), the formation of polycyclic 
aromatic hydrocarbon (PAH) dust in carbon-rich AGB stars, the H2 formation on 
dust grains, and the H2 photo-dissociation due to far ultra-violet (FUV) light in a 
self-consistent manner. We focus mainly on disk galaxies with the total masses rang- 
ing from lO^^M© and lO^^M© in this preliminary study. The principle results are as 
follows: The star formation histories of disk galaxies can be regulated by the time evo- 
lution of interstellar dust, mainly because the formation rates of H2 can be controlled 
by dust properties. The observed correlation between dust-to-gas-ratios {D) and gas- 
phase oxygen abundances {Aq = 12 + log (0/H)) can be reproduced reasonably well 
in the present models. The disks show negative radial gradients (i.e., larger in inner 
regions) of H2 fraction (/hs)? PAH-to-dust mass ratio (/pah), D^ and Aq and these 
gradients evolve with time. The surface-mass densities of dust (Hdust) are correlated 
more strongly with the total surface gas densities (Sgas) than with those of H2 (^Hs)- 
Local gaseous regions with higher D are more likely to have higher /h2 in individual 
disks and total H2 masses (Mhs) correlate well with total dust masses (Mdust)- More 
massive disk galaxies are more likely to have higher D, /pah, and /hs and smaller 
dust-to-stellar mass ratios (i^dust = ^dust/^star)- Early- type E/SO galaxies formed by 
major galaxy merging can have lower i^dust than isolated late- type disk galaxies. We 
also compare between galactic star formation histories in the metallicity-dependent 
and dust-dependent star formation models and find no major differences. Based on 
these results, we discuss the roles of dust in chemical and dynamical evolution of 
galaxies. 

Key words: ISM: dust, extinction - galaxies:ISM - galaxies:evolution - in- 
frared:galaxies - stars: formation 



1 INTRODUCTION 

One of the many important roles of interstellar dust in galax- 
ies is the formation of molecular hydrogen (H2) on its surface 
(e.g., Gould & Salpeter 1963; Hollenbach k Salpeter 1971; 
Cazaux & Tielens 2002). Giant molecular clouds (GMCs) 
composed of H2 are the major formation sites of stars in 
galaxies (e.g., Blitz 1993; Fukui & Kawamura 2010). Global 
star formation rates of galaxies are observed to be well cor- 
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related with surface mass densities of interstellar gas (e.g., 
Schmidt 1959; Kennicutt 1998). Dust can originate from 
stellar winds of AGB stars (e.g., Ferrarotti & Gail 2006; 
Zhukovska et al. 2008) and supernovae (e.g., Kozasa et al. 
1991; Nozawa et al. 2003), and therefore the production rate 
of dust in a galaxy can be determined by the star forma- 
tion history that controls the formation rates of AGB stars 
and supernovae. Thus, the evolution processes of dust, gas, 
and stars are mutually related, and detailed investigation of 
these coevolution processes can lead us to the better under- 
standing of galaxy formation and evolution. 
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Recent observational studies by infrared space tele- 
scopes have revealed physical properties of dust, their spatial 
distributions, and their correlations with their host galaxy 
properties in nearby and distant galaxies (e.g., Draine et al. 
2007; Meixner et al. 2010; Roman-Duval et al. 2010; Takagi 
et al. 2010; Dunne et al. 2011; Kaneda et al. 2011; Cortese 
et al. 2012; Skibba et al. 2012; Smith et al. 2012). For ex- 
ample, Draine et al. (2007) have investigated the total dust 
mass, the mass fraction of the dust contributed by PAHs, 
and the correlation between gas-phase oxygen abundance 
(Ao=12 +log(0/H)) and the dust-to-gas ratio (D) for 65 
nearby galaxies and found that the average PAH fraction in 
galaxies with Aq > 8.1 is 3.55%. Smith et al. (2012) have 
found a significant difference in the dust-to-stellar-mass ra- 
tio between late-type spirals and SOs and thus provided a 
new clue to the origin of SOs. Meixner et al. (2010) and 
Skibba et al. (2012) have derived the detailed 2D maps of 
D, dust temperature, and PAH fraction in the Large Mag- 
ellanic Cloud (LMC). 

Recent multi-wavelength observational studies of 
nearby galaxies have investigated physical correlations be- 
tween surface densities of star formation (Ssfr), neutral 
hydrogen (H I; Shi), and H2 (Shs) and discussed key pa- 
rameters that control global star formation in galaxies (e.g., 
Bigiel et al. 2008; Leroy et al. 2008). These observation have 
also provided data sets for radial gradients of H i and H2 
gas, which are quite useful for discussing what can deter- 
mine the H I and H2 gas mass fractions in galaxies. High- 
resolution observational studies of GMCs for galaxies in the 
Local Group (e.g., LMC and M33) have provided vital clues 
to how individual GMCs composed of H2 form from local 
H I and how star formation proceeds within GMCs (e.g., 
Rosolowsky et al. 2003; Kawamura et al. 2009). A number 
of recent observational studies have investigated correlations 
between properties of dust (e.g., D and Sdust) and gas (e.g., 
chemical abundances and Shs) in galaxies (e.g., Leroy et al. 
2011). 

These observational results have raised the following 
three key questions on the origin of dust and gas in galaxies. 
The first is what physical processes are responsible for the 
observed spatial variation of dust properties within galax- 
ies (e.g., radial variation of D and PAH-dust mass fraction; 
Meixner et al. 2010). The majority of previous theoretical 
studies of dust formation and evolution in galaxies are based 
on one-zone (or multi-zone) chemical evolution models (e.g., 
Dwek 1998, D98; Lisenfeld & Ferrara 1998; Hirashita 1999; 
Edmunds 2001; Calura et al. 2008; Piovan et al. 2011). Al- 
though such previous models provided theoretical explana- 
tions for a number of key observational results on dust prop- 
erties in galaxies (e.g., Aq — D relation), they did not allow 
astronomers to discuss the spatial distributions of dust prop- 
erties and their correlations with other galaxy properties 
(e.g., Hubble morphological types). Therefore, it is largely 
unclear in these previous studies what physical mechanisms 
are responsible for the observed spatial distributions of dust 
properties in galaxies. 

The second question pertains to the determinant fac- 
tors for H2 properties and their correlations with galaxy 
properties (e.g., luminosities and metallicities) . Previously, 
Elmegreen (1993) theoretically discussed the importance of 
the pressure and radiation field of the interstellar medium in 
the transition from H to H2 in galaxies. Blitz & Rosolowsky 



(2004) also investigated the importance of gas pressure in 
H — H2 transition based on observational data from nearby 
28 galaxies. Recent numerical simulations have incorporated 
the H2 formation process from H I (e.g., Pelupessy et al. 
2006, P06; Robertson k Kravtsov 2008, RK08; Gnedin et 
al. 2009; Christensen et al. 2012), which means that the- 
oretical studies will soon contribute to the solution of the 
above problem. The H2 formation model dependent only on 
gaseous metallicities and densities proposed by Krumholz, 
McKee & Tumhnson (2009, KMT09) is practically useful 
for investigating H2 formation and star formation from H2 
in galaxies and thus used in recent semi-analytic models 
(e.g., Fu et al. 2010; Lagos et al. 2012) and numerical simu- 
lations (Kuhlen et al. 2012, K12). However, these recent H2 
formation models assume that dust abundances are linearly 
proportional to metallicities (i.e., dust-to- metal ratio, Dz, 
is constant), which is inconsistent with observations (e.g., 
Galamex et al. 2011) which show a large dispersion in D 
for a given Aq (i.e., significantly different Dz). Furthermore, 
such an assumption of constant Dz is incompatible with dust 
evolution models that predict significant Dz evolution with 
time in galaxies (e.g., Inoue 2003; Calura et al. 2008). Thus 
we need to improve H2 formation models by incorporating 
dust formation and evolution in order to discuss the second 
question in a more quantitative manner. 

The third question is on the origin of the observed cor- 
relations between dust and galaxy properties (e.g., dust-to- 
stellar mass ratio along the Hubble morphological types; 
Cortese et al. 2012) and between dust and gas properties 
(e.g., gas surface densities dependent on dust surface densi- 
ties; Leroy et al. 2011). In order to discuss this third question 
in a quantitative manner, theoretical studies need to model 
both (i) the time evolution of dust/gas contents and star 
formation rates in galaxies and (ii) the dynamical evolution 
of galaxies in a self- consistent manner. Only a few previous 
models included the formation and destruction of dust and 
H2 formation on dust grains in a self- consistent manner and 
thereby discussed the time evolution of dust properties (e.g., 
Hirashita & Ferrara 2002; Yamasawa et al. 2011). Since they 
did not explicitly include dynamical evolution of galaxies in 
their models, the above-mentioned correlations were not ad- 
dressed at all. 

Chemodynamical simulations are useful and powerful 
tools for investigating both the spatial distributions of stel- 
lar and gaseous abundances and galaxy scaling relations 
simultaneously (e.g., Theis et al. 1992; Bekki & Shioya 
1999; Kawata 2001; Revaz & Jablonka 2012; RJ12) and 
thus should be also useful for theoretical studies on spa- 
tial distributions and scaling relations of dust and H2 prop- 
erties: we can not discuss the above key problems related 
to dust and H2 distributions without performing chemody- 
namical simulations. However, even the latest chemodynam- 
ical simulations with more sophisticated models for the time 
evolution of variously different chemical abundances (e.g., 
Rahimi & Kawata 2012, RK12; Bekki et al. 2012; B12) did 
not explicitly include the formation and evolution of dust. 
Although recent hydrodynamical simulations have incorpo- 
rated H2 formation on dust grains (e.g.,P06), they did not 
incorporate the time evolution of dust properties in a self- 
consistent manner. Thus it is high time for us to develop a 
new chemodynamical model with dust and H2 formation and 
thereby to try to answer the above three key questions. 
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Table 1. A summary for recent observations to be compared with the present simulations. 

Properties Physical meanings 

Aq — D Correlations between gas-phase oxygen abundance (Aq) and dust-to-gas ratios (D) 

^O ~ /pah Correlations between Aq and PAH-to-dust ratios (/pah) 

/h2 — D The dependences of /h2 on D in local ISM 

^dust — ^gas Correlations between local surface densities of dust (Sdust) cind those of total gas (Sgas) 

^dust ~ ^H2 Correlations between local surface densities of dust (Sdust) ^ind those of H2 (5]gas) 

dD/dR Radial gradients of D 

dfpAu/dR Radial gradients of /pah 

^dust — Mgas Correlations between total dust masses (Mdust) and total gas masses (Mgas) 

^dust ~ ^star Correlations between dust-to-star mass ratios (i^dust) cind total stellar masses (Mgtar) 

^dust — ^star The dependences of local i^dust on local stellar surface densities (Sstar) 



References ^ 
Galametz et al. (2011) 
Draine et al. (2007) 
Prediction 
Leroy et al. (2011) 
Leroy et al. (2011) 
Pappalardo et al. (2012) 
Meixner et al. (2010) 
Corbein et al. (2012) 
Cortese et al. (2012) 
Prediction 



^ 'Prediction' means that only predicted properties are presented: observational results are yet to be obtained. Only one representative 
reference is given. A full list of relevant papers is given in the main text. 



The three purposes of this paper are as follows. Firstly, 
we describe the details of the new simulation methods by 
which we can investigate 3D distributions of dust, H2, and 
stars of galaxies in a self-consistent manner. Secondly, we 
present the preliminary results on the time evolution of dust 
and H2 contents and star formation histories of disk galax- 
ies and their dependences on the model parameters of the 
simulations. Thirdly, we discuss recent observational results 
for dust and H2 properties of galaxies and their correlations 
with their host galaxy properties based on the present sim- 
ulation results; several of these observations to be compared 
with the present simulations are summarized in Table 1. It 
should be stressed that these observations of the spatial dis- 
tributions of dust and H2 and their correlations with galaxy 
properties can be appropriately addressed only by chemody- 
namical simulations with dust evolution such as the present 
one. This is the first paper of a series of papers on the co- 
evolution of dust, gas, and stars in galaxies. We therefore 
focus exclusively on the details of the new numerical meth- 
ods and the preliminary results. We discuss other important 
issues such as the evolution of dust composition, the im- 
portant influences of dust on star formation and chemical 
evolution histories of galaxies, and the derivation of spec- 
tra energy distributions from the results of chemodynamical 
simulations in our future papers. 

The plan of the paper is as follows: In the next sec- 
tion, we describe our new chemodynamical model with the 
formation and evolution of dust and H2. In §3, we present 
the numerical results on the long-term evolution of phys- 
ical properties of dust and H2 and their correlations (e.g., 
Aq — D relation) in disk galaxies. In this section, we also dis- 
cuss the dependences of the results on the adopted model 
parameters. In §4, we discuss the latest observational re- 
sults on dust and H2 properties of galaxies derived mainly 
from the Herschel and Spitzer telescopes. We also compare 
the predicted dust scaling relations with the corresponding 
observed relations in this section. We summarize our con- 
clusions in ^5. 



2 THE CHEMODYNAMICAL MODEL 

We mainly investigate the time evolution of dust and H2 
properties in forming disk galaxies embedded in massive 
dark matter halos with their physical properties (i.e., radial 



density profiles) consistent with predictions of a Cold Dark 
Matter (CDM) cosmology (e.g., Navarro et al. 1996; NFW). 
The present chemodynamical model incorporates chemical 
evolution of variously different elements (e.g.. He, C, N, O, 
Mg, and Ca), chemical enrichment by Type la. Type II, 
and aspherical supernovae (SNIa, SNII, and ASN respec- 
tively) and AGB stars, supernova feedback effects of SNIa 
and SNII, formation and destruction of dust, H2 formation, 
metallicity-dependent radiative cooling, and H2-regulated 
star formation. Therefore the present model is a much im- 
proved version of those used in investigating chemodynam- 
ical evolution of galaxy mergers (Bekki & Shioya 1999) and 
disk galaxy evolution in groups (Bekki & Couch 2011). 

In order to perform numerical simulations on GPU clus- 
ters, we have revised our previous chemodynamical code 
('GRAPE-SPH'; Bekki 2009) that can be run on the special 
computer for gravitational dynamics (GRavity PipE; Sugi- 
moto et al. 1990). In the present paper, we describe only the 
key ingredients of the code; the full details of the new code 
(including the code performance) will be described in our 
forthcoming paper (Bekki 2013). Recent numerical simula- 
tions with H2-regulated star formation (P06; R08; K12), the 
results of which are compared with those of this paper, do 
not include the important feedback effects of ASN on the in- 
terstellar mediums (ISM) of galaxies, which are investigated 
in B12. We therefore disable the code's function of feed- 
back effects and chemical enrichment by ASN in the present 
study. For convenience, physical meanings of symbols often 
used in the present study are summarized in Table 2. 



2.1 Chemical enrichment 

Chemical enrichment through star formation and metal ejec- 
tion from SNIa, II, and AGB stars is considered to proceed 
locally and inhomogeneously. SNe and AGB stars are the 
production sites of dust, and some metals ejected from these 
stars can be also accreted onto dust grains in the ISM of 
galaxies. We investigate the time evolution of the 11 chemi- 
cal elements of H, He, C, N, O, Fe, Mg, Ca, Si, S, and Ba in 
order to predict both chemical abundances and dust prop- 
erties in the present study. The mean metallicity Z for each 
A;th stellar particle is represented by Zk. The total mass of 
each jih (j = 1 — 11) chemical component ejected from each 
kth stellar particles at time t is given as 
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' ma,kYk,j{t — ifc), 



(1) 



where rris^k is the mass of the /cth stellar particle, Ykj{t — tk) 
is the mass of each jth chemical component ejected from 
stars per unit mass at time t, and tk represents the time 
when the A;th stellar particle is born from a gas particle. 
Azl^ At) is given equally to neighbor SPH gas particles (with 
the total number of Nnei,k) located around the A;th stellar 
particle. Therefore, the mass increase of each jth chemical 
component for ith gas particle at time t {Az^j{t)) is given 



^4j(^)= Y1 ^-,kYk,j{t-tk)/Nn, 



(2) 



where A^nei,* is the total number of neighbor stellar particles 
whose metals can be incorporated into the ith gas particle. 
We consider the time delay between the epoch of star 
formation and those of supernova explosions and commence- 
ment of AGB phases (i.e., non-instantaneous recycling of 
chemical elements). Therefore, the mass of each jth chemical 
component ejected from each ith stellar particle is strongly 
time-dependent. In order to derive the mass-dependent life- 
times of stars that become SNe la, SNe II, and AGB stars, we 
estimate the main-sequence turn-off mass (ttito) for stellar 
particles. We do so by using the following formula (Renzini 
& Buzzoni 1986): 



logmTo(ts) ^ 0.0558(logts)^ - 1.338 log ts + 7.764, 



(3) 



where mxo is in solar units and time ts in years. 

We adopt the 'prompt SN la' model in which the delay 
time distribution (DTD) of SNe la is consistent with recent 
observational results by extensive SN la surveys (e.g., Man- 
nucci et al. 2006; Sullivan et al. 2006). In this prompt SN 
la mode, there is a time delay (tia) between the star forma- 
tion and the metal ejection for SNe la. We here adopt the 
following DTD {g{ti^)) for 0.1 Gyr ^ tia ^ 10 Gyr, which 
is consistent with recent observational studies on the SN la 
rate in extra-galaxies (e.g., Totani et al. 2008; Maoz et al. 
2010, 2011): 



gia{tl; 



^g^T 



(4) 



where Cg is a normalization constant that is determined by 
the number of SN la per unit mass (which is controlled by 
the IMF and the binary fraction for intermediate- mass stars 
for the adopted power-law slope of —1). This adoption is 
pointed out to be necessary to explain the observed chemical 
properties of the LMC (Bekki & Tsujimoto 2012) 

The fraction of the stars that eventually produce SNe 
la for 3-8M0 has not been observationally determined and 
thus is regarded as a free parameter, /b. It is confirmed that 
the present results on dust and H2 properties do not depend 
so strongly on /b for 0.03 ^ /b ^ 0.09 in the present study. 
We therefore show the results of the models with /b = 0.09, 
which is a reasonable value for investigating luminous disk 
galaxies like the Galaxy (e.g., Tsujimoto et al. 2010). The 
chemical yields adopted in the present study are the same as 
those used in Bekki & Tsujimoto (2012) except those from 
AGB stars. We adopt the nucleosynthesis yields of SNe II 
and la from Tsujimoto et al. (1995; T95) and AGB stars 
from van den Hoek & Groenewegen (1997; VG97) in order 
to estimate Yk,j{t — tk) in the present study. 



Table 2. Description of physical meanings for symbols often 
used in the present study. 



Symbol 


Physical meaning 


Ao 


gas-phase oxygen abundances 


D 


dust-to-gas ratio 


fa 2 


mass fraction of molecular hydrogen (H2) 


/pah 


PAH-to-dust mass ratio 


^H2 


surface mass density of H2 


Sgas 


surface mass density of gas (H 1+H2) 


Sz 


surface mass density of metals (Z) 


^dust 


surface mass density of dust 


^PAH 


surface mass density of PAH dust 


^star 


surface mass density of stars 


MH2 


logio Sh2 


Mdust 


loglO ^dust 


Mg 


loglO ^gas 


/is 


loglO ^star 


^dust 


mass ratio of dust to star (M^ust/^star) 


-^dust 


initial dust-to-metal ratio in gas 


/dust,m 


mass fraction of metals locked up in dust 


Pth 


threshold gas density for star formation 


Tacc 


dust accretion timescale 


'Tdest 


dust destruction timescale 


Mh 


initial dark halo mass 



2.2 Dust model 

The present dust model is essentially the same as that 
adopted in the previous multi-zone model by D98, which 
reproduced reasonably well the observed chemical and dust 
properties of the Galaxy in a self-consistent manner. The 
dust model consists of the following four components: (i) 
production in stellar winds of SNe la and SNe II and AGB 
stars, (ii) accretion of metals of ISM on dust grains, (iii) 
destruction of dust by energetic SN explosions, and (iv) 
PAH formation. The present model is somewhat idealized 
in that it does not include coagulation of small dust grains 
and time evolution of dust sizes that have been recently in- 
vestigated in some one-zone models (e.g., Hirashita 2012). 
In the present paper, we focus exclusively on the most im- 
portant processes related to dust formation and evolution 
in galaxies. The influences of dust coagulation and dust size 
evolution will therefore be investigated in our forthcoming 
paper. 



2.2.1 Production 

The total mass of jth component (j=C, O, Mg, Si, S, Ca, 
and Fe) of dust from kih type of stars (k = I, II, and AGB 
for SNe la, SNe II, and AGB stars, respectively) is described 
as follows; 



k 
"^dustj' 



— ^^ F ■ 



/ k 

(me 



(5) 



where 6^^^ is the condensation efficiency (i.e., the mass frac- 
tion of metals that are locked up in dust grains) for each jth 
chemical component from kth stellar type, Fej is the func- 
tion that determines the total mass of metals that can be 
used for dust formation, and Tn^j j is the mass of jth compo- 
nent ejected from kth stellar type. The total mass of stellar 
eject a is estimated by using stellar yield tables by T95 and 
VG97. For stars with the initial masses (??ii) larger than 
8M0, m^dnst,,j is as follows: 
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Figure 1. The H2 mass fraction (/H2) of ^i gas particle with the SPH smoothing length (h) of lOOpc as a function of D (i.e., dust-to-gas 
ratio) for different x (ISRF strength) and n (hydrogen number density): x = 1 and n = 1 cm~^ (upper left), x = 1 and n = 10 cm~^ 
(upper right), x = 10 and n = 1 cm~^ (lower left), and x = 10 and n = 10 cm"*^ (lower right). The constant Do is the dust-to-gas 
ratio of the solar neighborhood in the Galaxy. It should be noted here that /h2 depends on h for a given x, n, and D, because hydrogen 
column densities (A/", Ni^ and A^2), which are key parameters for /h2: ^ire calculated from n (ni and 722) and h: the results shown in this 
figure are true only for h = lOOpc. Although h is fixed for different n and x in this figure, h is different for different n (h gc n-*^/^) and 
X in real simulations. Accordingly, this figure is used for an illustrative purpose (to demonstrate the important role of D in determining 



II 



16 E 



1 = 1 "' 



'5^ej,fcM 



for other than O 
for O 



(6) 



In the above estimation of m^ust,©? ^^e summation is done 
for / = Mg, Si, S, Ca, and Fe (i.e., ni — 5) and /j^i is the 
mass of Ith element in atomic mass units. This formula for 
SNe II is used for SNe la in the present study. Slj and Sl]j 
are 0.8 for j=Mg, Si, S, Ca, and Fe and 0.5 for c'. 

For stars with mi ^ SM© and C/O > 1 in the ejecta. 



''dust,, J 



AGB 

'^dust,j 



{f 



c«°c'-0.75m,^°J') forC 



for other than C 



(7) 



The ejecta from these C-rich stars is assumed to have 
^^c^ = 1 in the present study (i.e., the same value adopted 
by D98). 

For stars with mi ^ 8M0 and C/O ^ 1 in the ejecta. 



^^dusf,,j is as follows: 



AGB 



16 E 



ni cAGB II 



e],k/fJ^l 



for C 

for other than O & C (8) 

for O 



In the above estimation of m^S?,©? ^^e summation is done 
for / = Mg, Si, S, Ca, and Fe (i.e., m = 5). 6^f^ are 0.8 
for Mg, Si, S, Ca, and Fe. Although a range of S^j values 
could be adopted in the present study, we investigate mainly 
the models with the parameter values adopted above. These 
dust yield models are referred to as 'D98-type yield' for con- 
venience in the present study. We discuss the total amount 
of dust and do not discuss the dust composition in detail. 

For comparison, we also investigate some models in 
which the mass fraction (Fdust) of dust among all gas (met- 
als + dust) ejected from a star does not depend on the stellar 
mass. These models are referred to as 'uniform yield' for con- 
venience. We investigate models with uniform yield, mainly 
because several earlier one-zone models adopted such models 
for their investigation on the origin of the observed Aq — D 
relation. Although the dust mass fraction (Fdust) is a (fixed) 
free parameter, we show the models with Fdust =0.1. This is 
because the models with Fdust = 0.1 can better explain the 
observed D of galaxies. We mainly show the results of the 
models with D98-type yield rather than those with uniform 
one in the present study. 



6 K. Bekki 



Table 3. Description of the basic parameter values for the 
fiducial model. 



Physical properties 

Total Mass^ 

Structure _ 

Gas fraction 

Spin parameter 

Chemical yield 

Initial metallicity 

Dust formation 

PAH_^ 

Dust yield 

Initial dust/metal ratio 

Hubble-type 

Feedback £ 

SF_^ " 

IMF 

Particle number 

Softening length 

Gas mass resolution 



Parameter values 

r^ir = 120 kpc, c = 10 

/g = 0.1 

A = 0.038 

T95 for SN, VG97 for AGB 

[Fe/H]o = -3 

Tacc = 0.25 Gyr, Tdest = 0.5 Gyr 

i^PAH = 0.05 

D98-type 

0.1 

late-type disk (isolated) 

SNIa (/b = 0.09) and SNII 

H2-dependent, ISRF, pth = 1 cm~^ 

Salpeter {a = 2.35) 

A^dm = 900000, Ag = 100000 

^dm = 935 pc, Cg = 94 pc 

mg = lO^M© 



^ Mh = Mdm + Mg, where Mdm and Mg are the total masses 
of dark matter halo and gas in a galaxy, respectively. 

The NFW profile with a virial radius (rvir) and a c param- 
eter is adopted for the structure of dark matter halo. 
^ ^PAH is the mass fraction of PAH dust to total dust in the 
stellar eject a of C-rich AGB stars. 

^ /b is the binary fraction of stars that can finally explode as 
SNe la. 

^ Pth is the threshold gas density for star formation and in- 
terstellar radiation field (ISRF) is included in the estimation 
of H2 mass fraction in this model. 

2.2.2 Accretion 

Dust grains can grow by accretion of metals of ISM onto pre- 
existing cores and this accretion process is included in previ- 
ous models (D98). Following D98, we consider that the key 
parameter in dust accretion is the dust accretion timescale 
(Tacc). In the present study, this parameter can vary between 
different gas particles and is thus represented by Tacc,* for iih 
gas particle. The mass of jth component (j=C, O, Mg, Si, 
S, Ca, and Fe) of dust for iih gas particle at time t (dij(t)) 
can increase owing to dust accretion processes. The mass 
increase is described as 



Adtj(t) = AU{1 - fdust,^,J)d^,J{t)/7 



(9) 



where Ati is the individual time step width for the iih gas 
particle and fdust,i,j is the fraction of the jth chemical el- 
ement that is locked up in the dust. Owing to this dust 
growth, the mass of jth chemical component that is not 
locked up in the dust {zij{t)) can decrease, which is simply 
given as 



Azf^;{t) = -AU{l-U 



,)di,j{t)/Ts, 



(10) 



As is clear in these equations, the total mass of jth compo- 
nent in ith gas particle {mij{t)) is Zij{t) -\-dij{t). It should 
be stressed that the notation for dust mass here is for each 
particle whereas that for dust mass in §2.2.1 is for each star. 
Although Tacc,i can be locally different and time- 
dependent owing to complicated chemical and dynamical 
evolution of galaxies, we adopt a fixed value for all gaseous 
particles in a simulation. The fixed value (denoted Tacc) 



"^acc-O 13 Gyr ^ . 

T_ = 0.25 Gyr ^ 

^acc = 0-38 Gyr 
Tacc=0-50 Gyr 




2 + log(0/H) 



Figure 2. The time evolution of simulated galaxies on the 
Aq — D (upper) and Aq — /pah planes (lower) for four rep- 
resentative models with Tacc = 0.13 Gyr (solid, red), Tacc = 0.25 
Gyr (dotted, blue), Tacc = 0.38 Gyr (short-dashed, magenta), 
and Tacc = 0.50 Gyr (long-dashed, green). The gas-phase oxygen 
abundance (Aq = 12 + log (O/H)) is denoted as Aq in this figure. 
For these models, Mh = 10^^ M© and A = 0.038 are adopted. The 
large thick crosses indicate the observed location of the LMC on 
the two planes for comparison. The observational data are from 
Meixner et al. (2012) and Leroy et al. (2011), and the error bars 
of 0.2 dex are shown for Aq. 



varies between different models. As described later, we need 
to choose carefully this Tacc in order to reproduce the ob- 
served dust properties (e.g., D) in disk galaxies. Thus we 
avoid introducing additional free parameters for Tacc,i possi- 
bly dependent on time and the ISM properties in the present 
study. 



2.2.3 Destruction 

Dust grains can be destroyed though supernova blast waves 
in the ISM of galaxies (e.g., McKee 1989) and the destruc- 
tion process is parameterized by the destruction time scale 
('Tdest) in previous one-zone models (e.g., Lisenfeld & Fer- 
rara 1998; Hirashita 1999). Following the previous models, 
the decrease of the mass of jth component of dust for ith 
gas particle at time t due to dust destruction process is as 
follows 



AdfTit) = 



-Atidij{t)/Tde 



(11) 



where Tdest, i is the dust destruction timescale for ith particle. 
The dust destroyed by supernova explosions can be returned 
back to the ISM, and therefore the mass of jth chemical 
component that is not locked up in the dust increases as 
follows: 



AziJ^{t) = Atidij{t)/rdest,i 



(12) 
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Although rdest,i can possibly vary for different particles 
owing to different physical conditions of the local ISM, we 
adopt a fixed value (denoted as Tdust) for all gas particles in 
the present study. We consider that a reasonable Tdest is 0.5 
Gyr (D98) and show the results of the models with Tdest — 
0.5 Gyr: As described later (§2.7), only a narrow range of 
Tdest is allowed for reproducing observational results. The 
dust destruction can occur only for gas particles that are 
located in the surrounding regions of SNe la and SNe II 
only for a timescale of Tdest in the present study. 

Thus the equation for the time evolution of jih compo- 
nent of metals for iih gas particle are given as 

Zi,,{t + Ati) = Zij{t) + Azl^^it) + Azt:'-{t) + Azfj\t) (13) 

Likewise, the equation for dust evolution is given as 

di,,{t + AU) = d^,j{t) + AC;(i) + Adtfit) (14) 

Dust is locked up in stars as metals are done so, when gas 
particles are converted into new stars. This means that star 
formation process itself has an effect of destroying dust in 
the present study. 



2.24 PAH 

A growing number of observational studies on PAH proper- 
ties have been accumulated for galaxies within and beyond 
the Local Group (e.g, Draine et al. 2007; Meixner et al. 
2010; Takagi et al. 2010; Sandstrom et al. 2012). It is ac- 
cordingly timely for the present study to discuss the origin 
of the observed PAH properties in galaxies by using the new 
chemodynamical model. The most promising formation site 
of interstellar PAH dust is C-rich AGB stars, though di- 
rect observation supporting PAH formation in stellar winds 
of AGB stars is very weak (e.g., Tielens 2008 for a recent 
review) . We consider that some fraction of carbon dust pro- 
duced by C-rich AGB stars {C/O > 1) can finally become 
PAH dust and thereby investigate the PAH properties in the 
present study. The mass fraction of PAH dust to total car- 
bon dust in the eject a of C-rich AGB stars is a parameter 
denoted as Rpau • Using low-resolution simulations with dif- 
ferent i?pAH, we find that Rpau = 0.05 can better explain 
observations (described later in §2.7). We therefore show the 
results of the models with Rpau = 0.05 in the present study. 



2.3.1 H2 formation on dust grains 

The H2 formation on dust grains via grain catalysis has 
been extensively investigated by many authors (e.g., Gould 
& Salpeter 1963; Hollenbach & Salpeter 1971; Cazaux & 
Tielens 2002). A key parameter for the process is the total 
grain geometric cross section per H nucleon and is defined 
as 



idrigr 
da 



(15) 



where a, rigr, and n are the sizes of dust grains, the size 
distribution, and the hydrogen number density, respectively 
(Draine 2009). Sgr is often denoted in units of 10~^^ cm^ 
H"^ (Spitzer 1978) for convenience (i.e., Sgr = S_2i 10~^^ 
cm^ H~^). The H2 formation rate is given as 



dn2 
~dt 



— Rgrnm , 



(16) 



where ni, 712, and Rgr are the number density of atomic 
hydrogen, molecular hydrogen, and the effective H2 rate co- 
efficient, respectively. Rgr is given as 



l/8kBT 1/2/ . 

Rgr = -( ) (egr>Sg 



(17) 



2 vrmn 

where ks, T, mn, and egr are the Boltzmann constant, the 
temperature of gas, the mass of a hydrogen atom, and the 
formation efficiency averaged over the grain surface area, 
respectively. Observationally, Rgr is determined ('-^ 3x 10~^^ 
cm^ s~^; Jura 1975) so that the product of (cgr) and Sgr can 
be estimated ((egr)Sgr ~ 0.5; Draine 2009). 

We do not investigate the time evolution of the dust 
size distributions of galaxies in the present study. We as- 
sume therefore that (egr) is constant and thus Rgr is deter- 
mined by the dust-to-gas ratio (D). We adopt a reasonable 
value of (cgr) = 0.5 for classical silicate and carbonaceous 
grains (Draine 2009), though the adopted value could be a 
lowest possible value for the grains with the sizes larger than 
0.01/im (Draine 2009). Rgr for ith gas particle is determined 
as follows 



Rg- 



"-«""<iffi<i"''s'' 



(18) 



where Ti and Di are the temperature and the dust-to-gas 
ratio of ith gas particle, respectively, and Dq is the dust-to- 
gas ratio for the Galaxy (0.0064; Zubko et al. 2004). 



2.3 H2 formation and dissociation 

Previous different numerical simulations estimated the mass 
fractions of H2 to total hydrogen gas (/H2) in galaxies by 
using different H2 formation models (P06; R08; K12). The 
model adopted in the present study is essentially similar to 
that used in P06 in the sense that /h2 is determined by 
local far-UV (FUV) radiation fields and gas densities. One 
of the important differences between P06 and the present 
study is that the time evolution of dust abundances and 
compositions are explicitly followed and used for estimating 
H2 formation rates of local regions of galaxies in the present 
study. H2 for each local region in a galaxy is determined by 
the balance between H2 formation and H2 dissociation by 
FUV radiation. The formulas used in Goldshmidt & Stern- 
berg (1995) and Draine (2009) are adopted in deriving /h2- 



2.3.2 Formation/destruction equilibrium 

Goldshmidt & Sternberg (1995) adopted a plane-parallel 
equilibrium cloud model and thereby estimated the mass- 
ratios of H2 to H by solving the following equation that de- 
scribes the balance between H2 formation and dissociation 
by FUV radiation fields: 



RgrUm =XCfshield{N2)e ^712, 



(19) 



where C is the unattenuated H2 photo-dissociation rate for 
a unit incident FUV photon flux within the 11.2 to 13.6 
eV band, and x is the FUV intensity scaling factor relative 
to the unit FUV field, /shield (^^2) is the H2 self-shielding 
function, A^2 (A^i) is the column density of H2 (H i), and r 
is the optical depth for FUV continuum and given as r = 
aNi + 2crA^2, where a is the dust absorption cross section 
per hydrogen nucleus in the FUV wavelength range. 
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By substituting J nidr and J n2dr for Ni and A^2, 
respectively, in the above equation, a separate differential 
equation can be derived as follows 



RgrUdni = xC/shield(A''2)e 



(20) 



where n(=ni +712) is fixed. By solving the above equation, 
the total Hi column density (A^i°^) can be derived and writ- 
ten as a function of i?gr, a, Co, n, and C, and x as follows 

ATf*^ -\n{aaoGo + l), (21) 



where the dimensionless ao is described as 



and Go is described as 



Go^ [ 

Jo 



/shield (A^2)e 



-2(rN2 



dN2 



(22) 



(23) 



We can estimate x, ^gr, cr, and n for each gas particle by us- 
ing the local gas density, dust-to-gas ratio, and interstellar 
radiation field (ISRF) of the particle in the present simu- 
lations (described later). Therefore, by estimating Go for a 
reasonable function of /shield, we can derive A^i and N2 for 
each gas particle (by using N = Ni -\- 2N2 or n = ni + 2n2). 
Draine & Bertoldi (1996) derived an approximation for 
/shield as follows 

0.965 0.035 



/shield — 



(l+x/65)2 (1+X)0- 



■ exp[ 



1.5xl0-^(l+x)^-'],(24) 



where x = N2/5 x 10"*^ cm~^ and 65 = b/kms~^ (b is the 
Doppler broadening parameter and equal to 3 km s~^). Al- 
though this is a reasonably accurate approximation (Draine 
2009), we can not have an analytic formula for estimating 
Go (as a function of a), if we use this. Therefore, we adopt 
the following approximation from P06: 

/shield = (^)"°-'. (25) 

where A^ch is a characteristic column density of A^ch = 
1.75 X lO"*^"*^ cm~^. Since J^ e~^-y=dt — a/tt, we can derive 
an analytic formula for Go in equation (23), which allows 
us to save computational time significantly in the present 
simulations. 



2.3.3 ISRF 

In order to estimate x for each gas particle, we use the spec- 
tral energy distributions (SEDs) for stellar populations with 
different ages and metallicities derived by Bruzual & Chariot 
(2003; BC03). We first find stellar particles that are located 
within eg (i.e., gravitational softening length for gas parti- 
cles) of a gas particle, then derive the total flux at 1000 A 
around the gas particle from the sum of the SEDs of the stel- 
lar particles by using the SEDs of BC03. Finally, we estimate 
the FUV flux density around the gas particle for the local 
volume (VisRF oc /cfsRF^g? where /cisrf = 1) at each time 
step. We confirm that the present results do not depend on 
A^iSRF and discuss the results of the models with different 
A:iSRF in Appendix A. For consistency with the present IMF 
adopted for chemical evolution and star formation histories 
of galaxies, we use the SEDs of BC03 for the Salpeter IMF 
(a = 2.35). 



2.3.4 /h2 

We can estimate A^l,^ (and A^2,0 and thus /hs,* for ith gas 
particle by substituting the values of the temperature (T^), 
dust-to-gas ratio [Di)^ ISRF (xO? ^^^ hydrogen gas density 
(Ph) in these equations (19), (21), (22), (23), and (25). Since 
the H2 formation rate i?gr,i and dust absorption cross sec- 
tion {(Ti) depend on A (e.g., en — 1.9 x 10~^^ Di/Do cm^), 
/h2,2 depends strongly on dust properties in the present 
study. In converting Ni^i into ni,^ for ith gas particle, we 
use the local smoothing length {hi) of the particle (i.e., 
Ni^i — Jl ni^idr ~ 2r^l,^/^^). Fig. 1 shows an example of 
D-dependence of /h2 for a gas cloud size (h) of lOOpc, x = 1 
or 10, n = 1 or 10 cm"^ and five different T (=10, 30, 100, 
300, and 1000 K). 



2.3.5 Z - dependent YI2 and star formation 

Recent numerical simulations of galaxy formation and evo- 
lution with H2-regulated star formation (e.g., P06 and K12) 
assumed that dust abundances [D) of galaxies are linearly 
proportional to Z and the proportionality constant does not 
vary with time and location in galaxies. For convenience, the 
H2 and SF models adopted in the present study and that 
in P06 and K12 are referred to as D— dependent and Z- 
dependent, respectively, for clarity. An important question 
here is how different the predicted properties of galaxies are 
between numerical simulations with Z-dependent (e.g., P06 
& K12) and the present D-dependent H2 and star formation 
models. 

We investigate this question by adopting the same star 
formation recipe (i.e., Z-dependent) used in KMT09 for our 
disk galaxy formation model with Mh = 10^^ M© (the fidu- 
cial model). In estimating dust optical depth (tc) for local 
regions in a galaxy, K12 used the cell's column density in 
their adaptive mesh refinement simulation code. Since our 
code is not a mesh code, we estimate Tc for each local re- 
gion in a simulated galaxy by using the smoothing length 
of a SPH particle at the region. Since our main focus is not 
the Z-dependent model, we describe the results of only one 
Z-dependent model. 



2.4 Star formation 

We adopt the following two star formation (SF) recipes and 
compare between the results of models with the two recipes. 
One is the 'H-dependent' recipe in which star formation rate 
(SFR) around ith gas particle depends simply on the total 
gas density {pi). This has been a standard recipe since Katz 
(1992) implemented this in galaxy formation simulations. 
The other is 'H2-dependent', in which the SFR depends on 
H2 gas density rather than pi . This is more realistic and rea- 
sonable, given that star formation is ongoing within GMCs 
composed of H2 gas (e.g., K12). We mainly investigate mod- 
els with the H2-dependent SF recipe in the present study. 

In the H-dependent SF recipe, a gas particle is con- 
verted into a new star if (i) the local dynamical time scale 
is shorter than the sound crossing time scale (mimicking the 
Jeans instability) , (ii) the local velocity field is identified 
as being consistent with gravitationally collapsing (i.e., div 
v< 0), and (iii) the local density exceeds a threshold density 
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Figure 3. The time evolution of surface mass densities (in a logarithmic scale) projected onto the x-y (upper eight) and x-z (lower 
eight) planes for stars (/is, left) and gas (/ig, right). The time T is given in units of Gyr and a thick bar indicates 1 kpc. 



for star formation (pth)- Recent different models adopted dif- 

(RJ12), 1 cm-2 (RK12), 



(K12). We mainly inves- 
^, and briefly discuss 



ferent pth, for example, 0.1 cm 
10 cm"^ (R08), and 5 - 500 cm 
tigate the models with pth = 1 cm 
how the present results depend on pth later. Owing to the 
adopted SF model with pth, gas densities can not become 
very high. For example, the maximum gas density and min- 
imum dynamical time scale obtained for the fiducial model 
(described later) with a gas mass resolution of 10^ M© and 



Pth = 1 cm""^ are ^ 20 cm""^ and '-^ 3 x 10^ yr, respectively. 
The maximum density can be as large as ^ 8000 cm""^ in 
some models in which H2 formation is severely suppressed. 

Gas particles are converted into new stellar particles in 
the H2-dependent SF recipe described as follows. /h2 is esti- 
mated for each gas particle at each time step. A gas particle 
can be regarded as a 'SF candidate' gas particle if the above 
three SF conditions (i)-(iii) are satisfied. It could be possible 
to convert some fraction (oc /hs) of a SF candidate gas par- 
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Figure 4. The same as Fig. 3 but for H2 (AiH2 5 1^^) ^ii^d dust (/^dust^ right). 



tide into a new star at each time step until the mass of the 
gas particle becomes very small. However, this SF conver- 
sion method can increase dramatically the total number of 
stellar particles, which becomes numerically very costly. We 
therefore adopt the following SF conversion method. A SF 
candidate gas particle is regarded as having a SF probability 
of /h2- At each time step random numbers {R2; ^ R2 ^ 1) 
are generated and compared with /hs • If i?2 < /hs , then the 
gas particle can be converted into a new stellar one. 

Each iih stellar particle is born with a fixed IMF and 
an initial mass mo^i. The stellar mass decreases with time 



owing to mass loss 
the mass at time t 
mo, 2. The adopted 
Cimi~^, where mi 
and the slope a = 
The normalization 
mass cut-off), and 



a 



mo, 2 X (2 



rriu 



- ?7ir 



by SNe la, SNe II, and, AGB stars and 
(mi) can be significantly different from 
IMF in number is defined as '0(mi) = 

is the initial mass of each individual star 
2.35 corresponds to the Salpeter IMF. 

factor Ci is a function of mo,z, mi (lower 

TJiu (upper mass cut-off): 

^- (26) 



where m\ and niu are set to be O.IM© and lOOM©, respec- 
tively. We adopt a = 2.35 for all models in the present study. 
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Figure 5. The time evolution of star formation rates (SFRs, 
top), total stellar masses (Mgtar, middle), and gas mass fraction 



(Mgas/(Mg, 



•Msi 



bottom) for the two models with the 'H- 



dependent' (blue) and 'H2-dependent' (red) SF recipes. The red 
line corresponds to the fiducial model. The two models are exactly 
the same except the SF recipes. It should be noted that the gas 
mass fraction is estimated for all gas, including warm gas in the 
halo and well above the stellar disk: the gas mass fraction is fairly 
large (~ 0.6). The gas mass fraction is 0.38 in the stellar disk 
(R ^ 7.5 kpc and |2;| ^ 1 kpc), which is more consistent with the 
observed gas mass fraction of less luminous disk galaxies like the 
LMC (~ 30%; e.g., van den Bergh 2000). 



2.5 Gravitational dynamics, hydrodynamics, and 
feedback 

Since the present code is a revised version of our GRAPE- 
SPH code (Bekki 2009; Bekki & Couch 2011), the calculation 
of gravitational interaction between particles is based on the 
direct summation of gravitational force of the particles (i.e., 
no usage of tree codes). The detail and performance (on 
GPU clusters) of the new code will be given in our future 
papers (Bekki 2013). The direct gravitational calculation is 
done by the GPU whereas other calculations related to gas 
dynamic, star formation, and chemical evolution are done 
by the CPU in the GPU clusters used by the present study. 
One of key ingredients of the code is that the gravi- 
tational softening length (e) is chosen for each component 
in a galaxy (i.e., multiple gravitational softening lengths). 
Thus the gravitational softening length (e) is different be- 
tween dark matter (cdm) and gas (cg) and Cdm is determined 
by the initial mean separation of dark matter particles. Fur- 
thermore, when two different components interact gravita- 
tionally, the mean softening length for the two components 
is applied for the gravitational calculation. For example, 
e = (cdm + eg)/2 is used for gravitational interaction be- 
tween dark matter and gas particles. The total numbers for 



dark matter (A^dm) and gas particles (Ag) is 900000 and 
100000, respectively, in the fiducial model described later. 

In the present study, Cdm is relatively large (^ Ikpc 
for a model with Mh = 10^^ M©), which might severely 
suppress the formation of substructures in baryonic com- 
ponents. However, we think that owing to (i) the adopted 
multi-softening lengths and (ii) initially virialized dark mat- 
ter halo, such severe suppression of substructure formation, 
which could be regarded as a numerical artifact, is highly un- 
likely to occur. Indeed, clumpy structures of baryonic com- 
ponents can be formed in the present models, in particular, 
in H2 distributions. 

We consider that the ISM in galaxies can be modeled 
as an ideal gas with the ratio of specific heats (7) being 
5/3. The basic methods to implement SPH in the present 
study are essentially the same as those proposed by Hern- 
quist & Katz (1989). We adopt the predictor-corrector algo- 
rithm (that is accurate to second order in time and space) in 
order to integrate the equations describing the time evolu- 
tion of a system. Each particle is allocated an individual time 
step width (At) that is determined by physical properties of 
the particle. The maximum time step width (Atmax) is 0.01 
in simulation units, which means that Atmax = 1-41 x 10^ 
yr in the present study. Although a gas particle is allowed to 
have a minimum time step width of 1.41 x 10^ in the adopted 
individual timestep scheme, no particle actually has such a 
short time step width owing to conversion from gas to star 
in high-density gas regions. 

Each SN is assumed to eject the feedback energy (Esn) 
of 10^^ erg and 90% and 10% of ^sn are used for the increase 
of thermal energy ('thermal feedback') and random mo- 
tion ('kinematic feedback'), respectively. The energy-ratio 
of thermal to kinematic feedback is consistent with previous 
numerical simulations by Thornton et al. (1998) who in- 
vestigated the energy conversion processes of SNe in detail. 
The way to distribute ^sn of SNe among neighbor gas par- 
ticles is the same as described in B12. The radiative cooling 
processes are properly included by using the cooling curve 
by Rosen k Bregman (1995) for 100 ^ T < lO^K and the 
MAPPING III code for T ^ lO^K (Sutherland & Dopita 
1993). 



2.6 Initial conditions of galaxy formation 

Although we mainly investigate the evolution of forming disk 
galaxies embedded in massive dark mater halos, we also in- 
vestigate physical properties of merger remnants that are 
morphologically similar to early E/SO galaxies. This is be- 
cause we discuss the observed differences in dust and H2 
properties between different Hubble morphological types. 



2.6.1 Isolated disk galaxies 

Disk galaxies form from slow accretion of halo gas embedded 
in virialized massive dark matter halos in the present study. 
The halo gas of a forming disk galaxy is initially in hydro- 
static equilibrium determined by the mass distribution of the 
dark matter halo; the gas thereafter collapses and accretes 
onto the central region of the dark halo owing to radiative 
cooling processes. This disk formation model is essentially 
similar to those adopted by Kaufmann et al. (2007) for the 
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Table 4. A range of model parameters for the models with N = 10^. 



Parameters 
Range 



Mh (Mo)^ 
3 X 10^ - 10^ 



A_ Tacc (Gyr)^ SF recipe pth (atom cm ^) Dust yield _ 

0.02 -"o.06 0.13 - 0.5" H- or H2-dependent 1-10 D98-type or uniform 



'^ We mainly investigate the five representative models with M^ = 10^°, 3 x 10^°, 10^^,3 x 10^^, and lO^^M© for different A. We also 

investigate a model with M^ = 3 x lO^M© and A = 0.038 to search for a lowest possible D in low-mass disk galaxies. 

^ A = 0.02, 0.038, and 0.06 are investigated. 

^ Tdest is fixed at 0.5 Gyr for all models. 

^ D98-type yield is used for most models. For uniform dust yield model, F^ust (dust mass fraction among metals from SNe and AGB 

stars) is set to be 0.1. 




2 4 6 

Time (Gyr) 



Figure 6. The left two panes are for the time evolution of the forming disk galaxy in the fiducial model on the Aq — D and Aq — /dust,m 
planes, and the right four are for the time evolution of M^^ (top), M^ust (the second from top), fn^ (the second from bottom), and D 
(bottom). /dust,m is the mass fraction of metals (ejected from stars) that are locked up in dust. 



detailed investigation of angular momentum transfer in disk 
galaxy formation and by RJ12 for dwarf galaxy formation. 
The ratio of the dark matter halo mass (Mdm) to 
gaseous halo mass (Mg) in a forming disk galaxy is fixed 
at 9 for most models. The initial total halo mass is denoted 
as Mh (= Mdm + Mg). We adopt the NEW profile for the 
density distribution of the dark matter halo suggested from 
CDM simulations: 
Po 



p{r) 



(27) 



(r/rs)(l + r/rs)^' 

where r, po, and Vs are the spherical radius, the characteristic 
density of a dark halo, and the scale length of the halo, 
respectively. The c-parameter (rvir/^s, where rvir is the virial 
radius) is chosen for a given Mh according to the predicted 



c-Mh relation in the ACDM simulations (e.g., Neto et al. 
2007). 

The gaseous halo initially has the same spatial distri- 
bution as the dark matter and is assumed to be initially in 
hydrostatic equilibrium. The initial temperature of a halo 
gas particle is therefore determined by the gas density, total 
mass, and gravitational potential at the location of the par- 
ticle via Euler's equation for hydrostatic equilibrium (e.g., 
the equation lE-8 in Binney & Tremaine 1987). Therefore 
gaseous temperature Thaio{r) at radius r from the center of 
a disk galaxy can be described as: 



Thaio(r) : 



mp 1 / 

kB Phaio(r) J 



GM{r) 



■dr. 



(28) 
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Figure 7. The locations of gas particles on the Aq — D plane for four representative times steps in the fiducial model: T = lA Gyr 
(upper left), T = 5.6 Gyr (upper right), T = 9.9 Gyr (lower left), and T = 11.3 Gyr (lower right). The red square in each panel indicates 
the mean values of Aq and D. One of every 10 gas particles is shown in each panel. 



where rup G, and /cb are the proton mass, the gravitational 
constant, and the Boltzmann constant, respectively, and 
M{r) is the total mass within r determined by the adopted 
mass distributions of dark matter and baryonic components 
in the forming disk galaxy. 

The dark matter halo is dynamically supported by 
isotropic velocity dispersion and has no net angular mo- 
mentum in the present study. The gas has a net angular 
momentum with the radial distribution of specific angular 
momentum (j) described as j oc r that is consistent with the 
prediction of CDM simulations (e.g., Bullock et al. 2001). 

The spin parameter (A = ttttsTt) i^ a free parameter that 
controls the initial total amount of rotation and we mainly 
investigate A = 0.038 in the present study. The angular rota- 
tion of a gas particle (with the rotating axis being the 2;-axis) 
is determined by the adopted A, j(r), and the position of the 
particle. 



2.7 The choice of dust parameters 

Since no chemodynamical simulations have so far investi- 
gated variously different dust properties of galaxies, it is not 
so clear what a reasonable set of dust model parameters (i.e., 
Tacc and Tdest) is for explaining the observed dust properties. 
We accordingly investigate many 'low-resolution' {N ^ 10^) 



models with different parameter values of race and Tdest and 
thereby try to derive a reasonable set of the parameters that 
can explain observations. We then perform 'high-resolution' 
(N — 10^) simulations and investigate physical properties 
of gas, dust, and stars in galaxies by using the derived set of 
dust model parameters. The low-resolution simulations have 
N^m = 10^ and ATg = 2 X 10^. 

In determining reasonable Tacc and Tdest , we use the ob- 
served D and /pah of the LMC, because a growing number of 
observational data sets have been accumulated for the LMC 
(e.g., Meixner et al. 2010). We adopt the LMC-type models 
with Mh = 10^^ M0, A = 0.038 and variously different Tacc 
and Tdest and investigate the time evolution of D and /pah 
for 11.3 Gyr and compare the final D and /pah with the 
observations. Fig. 2 show the results for some of the models 
as well as the observed D and /pah with observational error 
bars. The observational data sets are from Meixner et al. 
(2010) and Leroy et al. (2011), and the error bars of 0.2 dex 
are shown for Ao, as done in Leroy et al. (2011) for possible 
spatial variation of Aq. In these models, Tdest = 0.5 Gyr and 
RvKR — 0.05. The models with Tacc = 0.13 Gyr and 0.25 Gyr 
can better explain the observed location of the LMG on the 
Ao — D plane than other ones. The model with Tacc = 0.13 
Gyr, however, shows /pah that is significantly lower than 
the observed value. The model with Tacc=0.25 Gyr shows 
/pah that is consistent with the observed value. 
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Figure 8. The radial gradients of fu^ (top), /pah (the second from top), D (the second from bottom), Aq (bottom) for the two 
representative time steps, T = 4.2 Gyr (left) and T = 9.9 Gyr (right) in the fiducial model. Small black dots indicate gas particles and 
red circles represent the mean values for each radial bin. One of every 10 gas particles is shown in each panel. 



Thus it is reasonable and realistic for the present study 
to adopt (racc,Tdest) = (0.25,0.5) Gyr in high-resolution sim- 
ulations. It is confirmed that as long as Tdest/Tacc ~ 2 (for 
Tdest = 0.13 — 0.5 Gyr), then the observed D and /pah can 
be reasonably well reproduced. Given that the above set of 
model parameters (race, T"dest) = (0.25, 0.5) Gyr) is consistent 
with those used for D98 that explains dust properties of the 
Galaxy self-consistently, we mainly investigate the models 
with (racc,T"dest) = (0.25,0.5) Gyr in the high-resolution sim- 
ulations. 



2.8 Parameter study 

We mainly describe the results of the 'fiducial' model for 
which the model parameters are described in Table 3. This 



fiducial model has a total mass of Mh = 10^^ M© that is sim- 
ilar to the mass of the LMC before its first passage of the 
Galaxy's virial radius (i.e., before the loss of its initial mass; 
Bekki 2011 and 2012). We adopt this fiducial model, because 
recent observational studies have provided a rich amount of 
information on dust and H2 properties of the LMC including 
2D distributions of D, /pah, and /hs across the LMC (e.g., 
Kawamura et al. 2009; Meixner et al. 2010), which can be 
compared with the fiducial model in detail. We also present 
the results of some representative models with different pa- 
rameters and the range of model parameters investigated in 
the present study are shown in Table 4. For all models, the 
initial iron abundances ([Fe/H]) and dust-to- metal ratios for 
gas particles are set to be —3 and 0.1, respectively. Owing 
to the adopted non-zero initial dust mass in gas, H2 for- 
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Figure 9. Correlations of dust surface densities (Sdust) with H2 surface densities (SH25 upper left), PAH surface densities (Spah, 
upper right), gas surface densities (Xlgas, lower left), and metal surface densities (Ez, lower left) at T = 11.3 Gyr (i.e., final time step) 
in the fiducial model. Each black dot represents the azimuthally averaged value at each radial bin. 



mation is possible from the very early stage of disk galaxy 
formation. 

We investigate the time evolution of dust and H2 prop- 
erties for the last '-^ 11 Gyr, which corresponds to the pe- 
riod of disk growth via gas accretion after virialization of 
the dark matter halos. The present simulations are there- 
fore somewhat idealized in the sense that they can not de- 
scribe the very early merging of sub-galactic clumps that 
formed stellar halos of disk galaxies. However, we consider 
that these simulations enable us to grasp some essential gra- 
dients of long-term dust and H2 evolution in galaxies. We 
will adopt more realistic initial conditions of galaxy forma- 
tion based on a ACDM cosmology and thereby investigate 
physical properties of dust and H2 in our forthcoming pa- 
pers. In the following, T in a simulation represents the time 
that has elapsed since the simulation started. 



3 RESULTS 

3.1 The fiducial model 

3.1.1 Time evolution of dust and H2 

Figs. 3 and 4 show the time evolution of 2D distributions 
of gas, new stars, H2, and dust (/ig, /is, MH25 and /Xdust, 
respectively) for the last 11.3 Gyr in the fiducial model. For 



clarity, the time evolution of surface mass densities (in a 
logarithmic scale) is shown for each component. As a gaseous 
disk forms through early rapid accretion of halo gas, new 
stars can form from the central high-density gaseous regions 
of the disk (T = 1.4 Gyr). Multiple explosions of SNe in the 
early burst phase can blow some fraction of the early gas 
disk so that numerous giant (kpc-sized) gaseous holes can 
be developed (T = 2.8 and 5.6 Gyr). As the disk grows by 
further gas accretion, star formation can occur in the outer 
part of the disk so that the disk size becomes larger (T = 5.6 
and 11.3 Gyr). The giant holes can be still conspicuous in 
the outer part of the final disk owing to the presence of 
SN explosions (T = 11.3 Gyr). The final stellar disk has 
a disk-like structure with a higher density surrounded by 
a diffuse halo- like or thick disk- like stellar component. A 
strong stellar bar can not be formed during the '-^ 11 Gyr 
evolution in this model. 

The formation of H2 is possible from high-density 
gaseous regions of the central part of the disk in the early 
phase of the disk formation {T — 1.4 Gyr). The simulated 
H2 gas shows clumpy structures, particularly, in later epochs 
(e.g., T — 2.8 and 5.6 Gyr), and a more compact distribu- 
tion in comparison with the total gas distribution. As stars 
form from H2 gas across the forming disk, a larger amount 
of dust can be produced so that the formation efficiency of 
H2 can increase. As a result of this, the H2 distribution be- 
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comes more widespread and a larger amount of H2 gas can 
be contained in the disk (T = 2.7 and 5.6 Gyr). Owing to 
the lack of a strong bar, the efficient gas-transfer to the cen- 
tral region of the disk through dynamical effects of the bar 
on gas does not occur. Consequently, a very strong concen- 
tration of H2 gas can not be seen in this model, though /h2 
is higher in inner regions of the disk. The simulated distri- 
bution of dust does not look so clumpy as that of H2 gas 
throughout the disk formation and evolution. Clearly, the 
gas disk has significantly smaller surface dust densities in 
the outer regions. 

As shown in Fig. 5, SFR can rapidly increase during 
the first 1 Gyr evolution to reach as high as IM© yr~^ in 
the disk. The SFR then decreases slowly to finally become 
less than O.IM© yr~^. In this lower mass disk model, SN 
feedback effects play a vital role in suppressing star for- 
mation within the disk. Comparison between the fiducial 
model and the one with H-dependent SF recipe suggests 
that there is no remarkable difference in the long-term SF 
histories between the two models. However, there are some 
significant differences in SF histories between low-mass mod- 
els (Mh ^ 10^° Mo) with H- and H2-dependent SF recipes, 
as described later in Appendix C. SFR can more violently 
change in the H-dependent SF model for the first few Gyr in 
comparison with the fiducial one. This violent change is due 
to the fact that stars can form in small high-density clumps 
in a bursty manner and then star formation can be truncated 
by strong supernova feedback effects. In the fiducial model 
with the H2-dependent SF, the SFRs in small high-density 
regions can not so steeply rise, because the H2 densities are 
not so high (in spite of the high total gas densities). Slower 
consumption of gas in the clumps therefore ends up with the 
less violent change of the total SFR in the fiducial model. 
There are no major differences in the final total stellar mass 
(Mstar) and gas mass fraction (Mg/(Mg + Ms)) between the 
models with H- and H2-dependent SF recipes. 

Fig. 6 shows the time evolution of the disk galaxy on 
the Ao — D and Aq — /dust,m planes and the time evolution 
of the total dust (Mdust) and H2 masses (Mhs) and mean 
D and /hs- Here the mean D, Aq, and /dust,m values aver- 
aged among all gas particles are shown. D of the disk evolves 
along a straight line {D oc Aq) until Aq ^ 8.6 and then very 
steeply increases to become D ^ 0.01. The origin of this 
steep rise in D is explained as follows. In the late evolution 
phase of the disk (12+logO/H > 8.6), the star formation 
rate can drop significantly owing to the lower gas density 
and thus the number of SN explosions becomes small. As a 
result of this, dust destruction by SNe can become much less 
efficient. On the other hand, the dust can still grow rapidly 
owing to accretion of gas-phase metals onto the pre-existing 
grains (i.e., at the expense of the metals). Accordingly the 
gas-phase abundance (Aq) can slightly decrease during this 
rapid D increase. On the other hand, /dust,m can very grad- 
ually increase until Aq ?^ 8.6 and then very steeply increases 
almost without changing Aq. 

Although H2 can be rapidly produced during the first 
0.5 Gyr evolution, the H2 gas can be efficiently consumed by 
the first burst of star formation in the forming disk. Mb.2 can 
slowly increase after the early high SFR phase (T < 2 Gyr) 
and takes a peak value (^ 8 x 10^ M©) around T = 6 Gyr. 
After the peak, both Mh2 and /h2 decrease slowly owing 
to gas consumption by star formation. The decrease for the 
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Figure 10. Upper panel: A correlation of the stellar surface 
densities (Sstar) with the dust-to-star mass ratios (i^dust = 
^dust/^star) in individual local regions at T = 11.3 Gyr in 
the fiducial model. Smaller black dots are for individual local 
regions and large red dots indicate the mean i^dust cit each radial 
bin. Lower panel: Plots of gas particles on the D — fn^ plane at 
T = 11.3 Gyr in the fiducial model. Red dots indicate the mean 
/h2 at each D bin. The red line connecting red dots means that 
H2 is higher for gas with higher D. One of every 10 gas particles 
is shown for clarity. About 81% of gas particles shown here have 
/h2 = (i.e., no molecular hydrogen). 



last 5 Gyr is more clearly seen in Mh2 in comparison with 
/h2- Both Mdust and D can slowly increase even after the 
peak formation phase of H2 owing to dust mass growth via 
accretion of gas-phase metals in ISM. 

Fig. 7 shows the distributions of gas particles on the 
Ao — D plane at four different time steps (T = 1.4, 5.6, 
9.9, and 11.3 Gyr). Only gas particles that are within the 
disk (z ^ 1 kpc) are plotted so that Aq — D relations for 
the disk can be investigated. The particles initially form a 
narrow line, the slope of which is significantly shallower than 
the observed slope of the Aq — D relation {D oc Aq; T = 1A 
Gyr). As time passes by, the Aq—D relation becomes steeper 
and gas particles become more widely spread on the Aq — D 
plane (T = 5.6 Gyr). A larger number of metal-poor parti- 
cles with Ao < 8 can have higher D (> 0.01), because they 
are located where SN explosions are rare and therefore dust 
can more efficiently grow by accretion of gas-phase metals. 
The majority of stars can be finally located within a nar- 
row line {D oc Aq) on the plane at T = 11.3 Gyr, though 
dispersions of D for a given Aq is large. 
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Figure 11. The time evolution of D (upper) and /h2 (lower) for 
six models with race =0.25 Gyr (blue), race =0.38 Gyr (red), Z- 
dependent SF recipe (green), H-dependent one (cyan), pth = 10 
cm~^ (magenta), and uniform dust yield (black). The blue line 
corresponds to the fiducial model. The basic parameters for initial 
conditions of galaxy formation (e.g., Mh and A) are exactly the 
same between the six models. 



3.1.2 Radial gradients 

Fig. 8 shows that the radial gradients of /hs, /pah, D, 
and Ao evolve significantly with time. Dust abundances can 
more rapidly increase in the inner disk regions owing to their 
more rapid star formation and chemical evolution histories, 
which end up with more rapid increases of H2 formation 
rates in the inner regions. Consequently, the radial gradi- 
ent of /h2 becomes steeper between T = 4.2 and 9.9 Gyr, 
and the H2 gradients are always negative (i.e., higher in 
the inner regions). A weak negative radial gradient of /pah 
(^ -6 X 10~^ kpc~^) can be seen at T = 4.2 Gyr, and 
this gradient becomes even weaker (^ — 2 x 10~^ kpc~^) at 
T = 9.9 Gyr. The /pah dispersion for a given R appears to 
be larger at T = 4.2 Gyr than at T = 9.9 Gyr. The almost 
lack of /pah gradient is discussed later in terms of the recent 
observational results of PAH dust in the LMC. 

The disk galaxy shows a negative radial D gradient at 
T=4.2 and 9.9 Gyr (- -2 x 10"^ kpc"^ and -4 x 10"'^ 
kpc~^, respectively) and the gradient becomes steeper at 
later epochs. A negative radial Aq gradient (i.e., gas-phase 
metallicity gradient) can be clearly seen at the two epochs; 
the gradient becomes shallower at later epochs. Dispersions 
of four properties, /hs, /pah, D, and Ao are large for a 
given radius, which reflects the fact that chemical evolution 
and star formation histories can be quite different between 
different local regions. The derived characteristic negative 
radial gradients of /h2, D, and Ao can be seen in other 
models, though the amplitudes of the gradients depend on 
model parameters. 



3.1.3 Dust-gas-star correlations 

Fig. 9 shows that there are clearly positive correlation of 
Sgas and Sh2 with Sdust in the final disk (T = 11.3 Gyr), 
though the Shs — Sdust correlation appears to be weaker. 
The almost constant Sgas/Sdust (^ 100) reflects the mean 



D of the disk and therefore the slope of the Sgas — 5]dust 
relation depends strongly on the dust abundances of sim- 
ulated disk galaxies (and thus on their chemical and dust 
evolution histories). Leroy et al. (2011) investigated correla- 
tions of surface densities of dust (Hdust) with those of total 
gas (Sgas) and H2 gas derived from CO for galaxies in the 
Local Group. They found that a stronger Sdust — Sgas cor- 
relation can be more clearly in the LMC in comparison with 
^^dust — Shs correlation. These observational results are con- 
sistent at least qualitatively with the present results shown 
in Fig. 9. 

Fig. 10 shows that there is no/little correlation between 
5]star (local stellar density) and i^dust = Mdust/Mstar (dust- 
to-star mass ratio) in local regions of the final disk: this could 
be regarded as a very weak correlation in the sense that local 
regions with higher Sstar have smaller Mdust/Mstar- This 
very weak correlation for log (E star) < 7.5 M© kpc~^ can be 
seen in other disk models. However, the remnants of gas- 
rich major mergers shows a slightly different Estar — i?dust 
relation for higher Estar (Bekki 2013), which could be an 
important difference between late- and early- types galaxies. 

Fig. 10 also shows that there is a positive correlation 
between /h2 and D in local regions of the final disk, though 
dispersion in /h2 is large for a given D. The derived larger 
/h2 for higher D (or higher D in more H2-rich ISM) is rea- 
sonable, given that a larger amount of H2 can be produced 
owing to higher Rgr for higher D in the present H2 formation 
model. This result has an important implication on the evo- 
lution of radial D gradients of galaxies in clusters of galaxies, 
which is discussed later. 



3.2 Parameter dependences 

Here we briefly summarize the governing dependences of 
dust and H2 properties, and their correlations with other 
key model parameters (e.g., Tacc and Mh). 



3. 2. 1 Dust parameters 

The time evolution and final values of D and /h2 depend 
strongly on Tacc- As shown in Fig. 11, although the very 
early evolution of D and /h2 (T < 1 Gyr) is not different 
between the two models with Tacc = 0.25 and 0.38 Gyr, the 
model with shorter Tacc (i.e., more rapid dust growth) shows 
more efficient production of dust and H2 and thus has larger 
values of D and /h2 • For these two models, a factor of ^ 0.5 
difference in Tacc ends up with a factor of ^ 3 difference 
in D and /h2- This strong Tacc-dependence of D and /h2 
implies that the modeling of dust growth is very important 
for the formation and evolution of dust and H2 : Tacc needs to 
be carefully chosen when /h2 in galaxies is investigated by 
using numerical simulations with dust. D98 suggested that 
ISM-phase averaged dust accretion timescale can be as long 
a ^ (1 — 2) X 10^ yr, which is roughly consistent with the 
adopted reasonable value of Tacc • 

Fig. 11 shows that there are no significant differences 
in the evolution of D and /h2 between the fiducial mod- 
els with D98-type and the model in which uniform dust 
yields are adopted (i.e., no PAH dust production), as long 
as Fdust = 0.1 is adopted as in previous one-zone models. 
Fig. 12 also shows that the final disks in the two models 
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Figure 12. The same as Fig. 8 but for different models with Z-dependent SF recipe (upper left), H-dependent one (upper right), 
Pth = 10 cm""^ (lower left), and uniform dust yield (lower right). In the H-dependent model, fn^ is not estimated at all (it is set to be 
0). The PAH dust evolution is not included in the uniform dust model (i.e., /pah = 0)- Thus 'inapplicable' is shown for the relevant 
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have very similar radial profiles of H2, D, and Aq (See Fig. 
8 for comparison). These results imply that the evolution 
of H2 and dust is not sensitive to the inclusion of of de- 
pendence to dust production rate and dust composition on 
stellar masses. A reasonable value of Fdust, however, should 
be carefully chosen in a simpler model with uniform dust 
yields. 
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3.2.2 SF recipes 

The models with Mh = IO^^Mq yet different SF recipes and 
SF threshold gas densities pth are investigated. The three 
principle results on the dependences of the time evolution 
and final values of D and /h2 on SF recipes, which are shown 
in Figs. 11 and 12, are summarized as follows. First, there 
are no major differences in the D evolution between the 
models with the D-dependent (i.e., fiducial model) and H- 
dependent SF recipes, which suggests that D evolution is 
relatively insensitive to the modeling of global star formation 
in galaxies. 

Second, although D evolution is essentially the same be- 
tween the models with the D-dependent (i.e., fiducial model) 
and Z-dependent SF recipes, the /h2 evolution is signifi- 
cantly different between the two in the sense that /h2 is 
significantly lower in the model with the Z-dependent SF 
recipe. The origin of the lower /h2 is discussed later in §4.1. 
Third, a factor of 10 difference in pth (1 ^ pth ^10 atom 
cm~^) can yield a factor of '^ 2 difference in fa^ in the 
present study. However, the evolution and the final value 
of D do not depend so strongly on pth' D is only slightly 
smaller in the model with higher pth owing to stronger sup- 
pression of star formation that produces dust. 
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Figure 13. The dependences of D (top), /pah (middle), and 
^dust (= Mdust/Mstar; bottom) on Mh for A = 0.038. 



3.2.3 Galaxy mass and spin parameter 

Fig. 13 shows that D is larger in more massive galaxies 
(i.e., larger Mh) for Mh ^ lO^^M©. This Mi, - D rela- 
tion is not clearly seen in more luminous disk galaxies with 
Mh > lO^^M©: the Mh - D relation is quite flat. The final 
/pah depends weakly on Mh such that /pah is larger for 
larger Mh . The derived Mh — /pah relation may well be ap- 
proximated as /pah oc Mh"^^. Since more massive galaxies 
have higher chemical abundances in the present study, this 
result means that galaxies with higher Aq can have higher 
/pah. High-mass galaxies have higher metallicities so that 
/h2 can be higher. The time evolution of /h2 for high-mass 
(Mh = lO^^Mo) and low-mass (Mh = lO^^M©) disk galaxies 
is given and briefly discussed in Appendix C. 

The dust-to-star mass ratio (i?dust) depends on Mh such 
that it is smaller for larger Mh. Also, the Mh — i^dust rela- 
tion appears to be different below and above Mh = 10^^ M©: 
^dust depends more strongly on Mh for Mh > lO^^M©. The 
final disks in the models with higher A (0.06) show lower 
Mh2, Mdust, /h2 5 and D, mainly because larger disks with 
lower mass densities can be formed in the models. Negative 
radial gradients of dust and H2 properties and dust-H2 cor- 
relations clearly seen in the fiducial model can be also seen 
in the models with different Mh and A. 

Fig. 14 shows that irrespective of galaxy halo masses 
(Mh), final Sgas positively correlate with Mdust (i-c, Sgas oc 
5]dust) in disk galaxies. The slopes of the correlations are 



steeper in galaxies with lower Mh owing to their lower dust- 
to-gas ratios, and the maximum Sdust is larger for larger Mh. 
Positive correlations between Sdust and Sh2 can be seen, 
though they are significantly weaker in comparison with the 
Sdust — 5]gas relations. The derived stronger correlations in 
the Sdust — Sgas relations are consistent with observations by 
Leroy et al. (2011). The comparative model with race = 0.38 
Gyr, in which dust growth is slower, shows a steeper slope 
in the Sdust — 5]gas relation and a smaller maximum value of 
5]dust. These results for the comparative model imply that 
dust accretion processes can be a key parameter that con- 
trols Sdust — Sgas and Sdust — SH2 relations in disk galaxies. 



4 DISCUSSION 

4.1 Necessary to model dust evolution in 
predicting star formation histories ? 

Previous chemical evolution models clearly showed that 
dust-to-metal ratios (Dz) of galaxies evolve with time espe- 
cially in the early evolution phases of galaxies (e.g., Hirashita 
1999; Inoue 2003; Calura et al. 2008). These previous models 
are one-zone, in which star formation models are rather ide- 
alized (e.g., non-inclusion of H2 formation models) and dy- 
namical evolution controlling star formation processes is not 
explicitly included. Therefore it is important for the present 
study to investigate how the mass fraction of metals that 
are locked up in dust (/dust,m) can evolve with time in the 
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present chemodynamical model. Fig. 14, showing the /dust,m 
evolution for four representative models with different Tacc 
and Mh, confirms that /dust,m can rapidly change in the 
early evolution of galaxies (T < 4 Gyr) , though the means 
by which /dust,m changes depends on the model parameters. 
This result strongly suggests that the total amount of dust 
in galaxies can not be linearly proportional to their metal- 
licit ies (Z) especially in early galaxy formation phases. 

Recent numerical simulations of galaxy formation and 
evolution with H2-regulated star formation (e.g., P06 and 
K12) assumed that dust abundances (D) of galaxies are lin- 
early proportional to Z and the proportionality constant 
does not vary with time and location in galaxies. The re- 
sults shown in Fig. 14 suggests that this assumption is not 
realistic, especially, in the early formation phases of galax- 
ies. An important question here is how different the pre- 
dicted properties of galaxies are between numerical simula- 
tions with Z-dependent (e.g., P06 & K12) and the present 
D-dependent star formation models. Fig. 15 shows the time 
evolution of star formation rates and dust properties (D and 
/pah) in the two models with Mh = 10^^ M© and D- and 



Z— dependent star formation recipes. Although the SFR can 
more violently change with time in the first ^ 2 Gyr evo- 
lution for the model with the Z-dependent SF recipe, the 
overall SFR is very similar between these two models. As a 
result of this, the time evolution of dust properties is very 
similar between the two models. This similarity in the time 
evolution of SFR, D, and /pah can be seen in models with 
greater Mh (^ 10^^ M©). 



These results imply that the Z-dependent SF recipe is 
reasonable and realistic enough to investigate the time evo- 
lution of SFRs (regulated by H2 formation and evolution) in 
galaxies, given that additional new parameters such as Tacc 
and Tdest should be considered in the D-dependent SF recipe. 
One might not have to adopt this recipe, unless one wishes to 
investigate dust properties such as radial gradients of D and 
/pah in detail. One of advantages of the present chemody- 
namical model is that we can discuss the observed gas-phase 
metallicities, dust properties, and dynamical properties of 
galaxies in a fully self- consistent manner. Thus we can make 
the most of the present chemodynamical model when we in- 
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vestigate correlations between physical properties of dust, 
gas, and stars in galaxies. 

However, as shown in Fig. 11, the final /h2 is signifi- 
cantly lower in the Z-dependent SF model than in the D- 
dependent one. The derived low /h2 ('^ 0.05) is not consis- 
tent with /h2 fractions observed in less luminous disk galax- 
ies like the LMC (e.g., /hs ~ 0.09 - 0.16; van den Bergh 
2000). Dust-to- metal ratios can increase during chemical 
evolution of galaxies so that H2 production on dust grains 
can become more efficient. Since the Z-dependent model 
does not include this evolutionary effect, it can underpredict 
/h2 in the present SPH simulations. This does not mean that 
the Z-dependent model has a problem in predicting /h2 • The 
time evolution of /h2 depends on the adopted Tc (i.e., op- 
tical depth) and the fixed initial dust-to-metal ratio in the 
Z-dependent SF model (KMT09). If different values of these 
two parameters are chosen in the Z-dependent model, the 
final /h2 niight well become more similar to that derived 
in the D-dependent model and thus to the observed one. 
Accordingly, we do not regard the significant /h2 difference 
between the two models as a serious problem: depending on 
which model we adopt, we need to carefully chose the model 
parameters (e.g., Tdust and Tc etc). 



4.2 Origin of dust scaling relations 

We here discuss the observed three key dust scaling rela- 
tions: correlations between Aq and D (e.g., Galametz et 
al. 2011; Leroy et al. 2011), between Mstar and i^dust = 
Mdust/Mstar (e.g., Corbelh et al. 2012; Cortese et al. 2012), 
and between Mh2 and Mdust (Corbelli et al. 2012). 



4.2.1 Ac 



D relation 



Recent observational studies have confirmed that galaxies 
with higher Aq are more likely to show higher D, which 
is approximated roughly as D oc Aq (e.g., Galametz et al. 
2011; Leroy et al. 2011). Although the observed Aq - D 
relation is not strong and the observational estimation of 
dust masses of galaxies has some uncertainties, it is impor- 
tant for the present study to confirm whether the relation 
can be well reproduced by the present new model with dust 
formation and evolution. Fig. 17 shows the locations of six 
simulated disk galaxies different Mh (= 3 x 10^, 10^°, 3 x 
10^°, 10^\3 X 10^\ lO^^Mo) on the Aq - D plsme. For com- 
parison, the location of the Galaxy in the one-zone model 
by D98 and those of dwarf galaxies at Aq — 8.2, 8.4, and 
8.6 in the best one-zone model with the x parameter being 
30 by Lisenfeld & Ferrara (1998) are plotted in this figure. 
Clearly, the simulated galaxies have a Aq — D relation 
very similar to the observed relation, though the simulated 
D is slightly smaller than the observed one for a given Aq- 
Given the possible observational errors in the estimation of 
D, this similarity suggests that the present new model can 
explain very well the observed Aq — D relation. The present 
simulations predict systematically higher D for a given Aq 
for lower Aq (< 8.4) in comparison with the one-zone mod- 
els. The observed large dispersion in D at lower Aq < 8.4 
for a given Aq has not been well reproduced by the present 
models and thus needs to be investigated in our future stud- 
ies. 
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Figure 15. The time evolution of /dust,m for different models 
with Mh = lO^^Mo and race = 0.25 Gyr (solid, blue), Mh = 
lO^^Mo and race = 0.38 Gyr (dotted, green), Mh = IO^OM© 
and race = 0.25 Gyr (short-dashed, red), and Mh = IQ-'^^M© and 
race = 0.25 Gyr (long-dashed, cyan). 



In the present model, supernova feedback effects play a 
vital role in determining global star formation histories of 
galaxies and thus the total amount of gas that is converted 
into new stars in galaxies. Supernova feedback effects can 
less strongly suppress star formation in more massive galax- 
ies owing to deeper gravitational potential wells. As a re- 
sult of this chemical enrichment, dust production can pro- 
ceed more efficiently in more massive galaxies. Therefore, 
more massive galaxies can finally have higher metallicities, 
D, and mass-ratios of final stars to initial gas. The simulated 
Ao — D relation due to this effectiveness of supernova feed- 
back is strongly dependent on Mh. Thus the present study 
is the first chemodynamical simulation that has successfully 
reproduced the observed Aq — D relation. 



12.2 Rdu 



Mstar relation 



Recently Cortese et al. (2012) have found an interesting dust 
scaling relation that galaxies with larger Mstar are more 
likely to have higher i^dust • They have also found that i^dust 
is systematically lower in early-type E/SO galaxies in com- 
parison with late- type disk galaxies. Although these corre- 
lations of i?dust with Mstar and galactic morphologies show 
larger dispersions, they could have valuable physical mean- 
ings of galaxy formation and evolution. We have thus inves- 
tigated the i?dust — Mstar relation for a number of represen- 
tative models with different parameters. In order to discuss 
the dust properties of early- type galaxies, we have taken the 
^dust — Mstar relation of early- type galaxies from our another 
work on E/SO formation by gas-rich mergers (Bekki 2013). 
Fig. 18 shows that (i) i^dust is lower for larger Mstar and (ii) 
this relation appears to be steeper for log (Mstar /M©) larger 
than 9.6. These simulated trends are qualitatively similar to 
the observed ones (Gortese et al. 2012), which implies that 
the present model is quite good at grasping some essential 
ingredients of the formation process of this i^dust — Mstar 
relation. However, the present disk formation models can 
not reproduce the observed galaxies with very low i?dust 
(logi^dust < —3), in particular, those with luminous galax- 
ies with log(Mstar/M0) > 10. 

This failure of the disk formation models may reflect 
the limitation of the models in which violent merger events 
are not explicitly included. However, the major merger mod- 
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Figure 16. The same as Fig. 5 but for the models with with 'D- 
dependent' (blue) and 'Z-dependent' (red) SF recipes. The blue 
lines correspond to the fiducial model. In the two models, the 
SFR at each local region is estimated from local H2 densities. 
However, the way to estimate H2 mass fraction at each local re- 
gion is different between the two. The Z-dependent SF recipe is 
described in the main text. 
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Figure 17. The locations of six models with different M^ yet 
same A (=0.038) from the present simulation (blue circles) on the 
Aq — D plane. For comparison, the results of one-zone models 
(D98; Lisenfeld & Ferrara 1998) are shown by red circles. The 
dotted line is the observed relation by Leroy et al. (2011). 



els, in which two late- type disk galaxies can be transformed 
into early-type E/SO ones, can show significantly lower i^dust 
(Bekki 2013). Fig. 18 shows i^dust of the remnants of three 
luminous major merger models with prograde-prograde or- 
bital configurations for Mh ^ lO^^M©. As shown in Fig. 18, 
the remnant with Mh = lO^^M© shows logi?dust < —3, 
which implies that early-type galaxies should have lower 
i^dust, if they are formed from major merging. It should be 
here noted that recent simulations of gas-rich major galaxy 
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Figure 18. The location of five late- type disk galaxies with 
different Mh (blue) and three luminous (Mh ^ IQ^^Mq) early- 
type galaxies formed from major merging (red) on the Mstar — 
i^dust plane. The results of major mergers are taken from Bekki 
(2013). 



mergers by Hayward et al. (2011) suggested that the low 
^dust of a galaxy is due simply to consumption of metal- 
enriched gas. A minor fraction of galaxies have extremely 
low RdMst (logi^dust < —4), which can not be reproduced by 
any model in the present study. These galaxies might have 
experienced some environmental effects (e.g., dust stripping 
by hot intra-cluster medium) or evaporation of dust by some 
hot radiation sources. The origin of these galaxies will be ex- 
plored in our forthcoming papers. 



12.3 Mh, 



Md 



relation 



Corbelli et al. (2012) have recently investigated the total 
masses of dust, H2, and Hi in late- type Virgo cluster galaxies 
and found Mdust — Mu, and Mdust — Mgas relations (where 
Mgas = Mhi + Mhs) for the galaxies. The best-fit relation 
between Mu, and Mdust and between Mgas and Mdust are 
described as 

log Mh2 = 0.77(±0.12) X log Mdust + 3.2(±0.9) (29) 

and 



log Mgas = 0.75(±0.09) X log Mdust + 3.9(=b0.6), 



(30) 



respectively. The metallicity-dependent CO-to-H2 conver- 
sion factor is used for the total H2 mass estimation of 35 
galaxies in deriving the above relations. Fig. 19 shows these 
two observational relations as well as the results from the 
present models with different Mh and A for comparison. 

Clearly, the simulated Mdust — Mgas relation is quite 
similar to the observed relation, although the locations of 
the simulated massive disk galaxies (Mdust ^ 8 x 10^ M© 
or Mh ^ 3 X 10^^ M0) are slightly above the observed re- 
lation. The locations of the simulated less massive galaxies 
with Mdust < 4 X 10^ M© on the Mdust — Mu, plane are ap- 
preciably (^0.5 dex) below the observed relation. However, 
given the observed large error (0.9 dex) for the Mdust — Mu, 
relation, the discrepancy between the simulation and the 
observation can not be serious. Instead, the present model 
appears to do a good job in reproducing the Mdust — Mu, 
relation as well as Mdust — Mgas- The present results also 
suggest that the Mdust — Mu, relation can become steeper 
for Mdust ^ lO'^M©. The slightly under-abundant Mu, for 
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Figure 19. Big red dots show the locations of 15 simulated disk galaxies with different Mh and A (=0.02, 0.038, and 0.06) on the 
-^dust ~-^H2 (upper) and M^ust — ^Hl+H2 planes (lower). For comparison, observational results of 35 galaxies from Corbelli et al. (2012) 
are shown by small black dots. The dotted lines are the observationally derived relations from Corbelli et al. (2012) and observational 
error bars (~ 0.1 dex) are shown by a large cross. 



a given Mdust derived in the present study would need to be 
re-investigated in our future simulations with more realistic 
initial conditions of galaxy formation. 

4.3 Radial gradients of D and /h2 

Recent observational studies have revealed radial gradients 
of AjjY (UV attenuation), dust surface densities, and D in 
galaxies (e.g., Boissier et al. 2004; Pappalardo et al. 2012; 
PI 2). Recent theoretical works have investigated the time 
evolution of radial gradients of D and Dz for different mod- 
els of dust growth and destruction (e.g., Mattsson et al. 
2012). The present study has clearly shown that luminous 
disk galaxies can have negative D gradients (i.e., higher D 
in inner regions), which is qualitatively consistent with the 
observational results by P12 for non- HI- deficient galaxies in 
the Virgo cluster of galaxies. The slightly steeper gradients 
of H2 in comparison with D in P12 appears to be also con- 
sistent with our model predictions. 

One of interesting recent results on radial D gradi- 
ents (PI 2) is that the D gradients appear to be differ- 
ent between galaxies with different Hl-deficiency parameters 
(defui)'- The parameter defui is defined as the difference 
between the observed H I mass (in logarithmic units) and 
that expected for a galaxy of the same linear diameter and 



morphological type in a comparison sample of isolated ob- 
jects (Roberts & Haynes 1994). Galaxies with defui < 0.4 
(i.e., non- deficient) show negative D gradients whereas those 
with defui > 0.7 (i.e., strongly Hl-deficient) show slightly 
positive ones (Fig. 13 in P12). Given that no model in the 
present study shows such positive D gradients, this observa- 
tion needs to be discussed in the context of environmental 
effects on galaxies. If H i gas in the outer parts of galax- 
ies can be more efficiently stripped by some environmental 
effects (e.g., ram pressure or tidal stripping) in comparison 
with dust, then D might well can increase in the outer parts 
after this stripping process. On the other hand, if the inner 
parts remain intact during the stripping process, D of the 
inner parts would not change. Therefore, this more efficient 
stripping of H i gas in the outer parts of galaxies could end 
up with positive radial D gradients. 

It appears to be theoretically unclear whether and why 
this more efficient stripping of H i gas could be possible. 
However, Fig. 10, which demonstrates a larger amount of 
dust in gas with higher H2 mass fractions, suggests the fol- 
lowing scenario for the more efficient H i stripping. A larger 
amount of dust in the outer parts of disks can be locked 
up in GMCs in comparison with H i gas (as shown in Fig. 
10). Therefore, if gas stripping process is more efficient for 
H I gas, then D could increase after gas stripping. Fig. 20 
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Figure 20. The cumulative mass distributions of total gas (solid, red), dust (dotted, blue), H I (short-dashed, magenta), and H2 
(long-dashed, green) at T = 11.3 Gyr in the fiducial model. 



shows the cumulative mass distributions for H i, H2, dust, 
and total gas (i.e., H 1+H2) in the fiducial model. The more 
compact H2 distribution in Fig. 20 strongly suggests that 
H I can be more efficiently stripped by tidal or ram pressure 
stripping in the outer parts of disks. Therefore, the above 
scenario seems to be reasonable and realistic. Given that no 
previous studies have so far investigated the above scenario, 
it is our intention in future studies to confirm that the origin 
of the observed positive D gradients is really related to more 
efficient H i stripping of disk galaxies in clusters of galaxies. 



4.4 Origin of PAH dust properties 

The present study has shown that the PAH-to-dust mass ra- 
tio (/pah) is higher for luminous (or massive) disk galaxies 
with higher metallicities (^o)- This correlation is qualita- 
tively consistent with recent observational results by Draine 
et al. (2007), which shows /pah (which is equivalent to their 
^pah) is 1% for metal-poor galaxies with Aq < 8.1 and 
3.55% for metal-rich galaxies with Aq ^8.1. Interestingly, 
their results show a lack of galaxies with higher /pah > 0.02 
and Ao < 8.1. This is also consistent with the present PAH 
formation models which show a rapid rise of /pah only after 
Ao > 8.1 (i.e., very steep slope of the simulated Aq — /pah 
relations for Aq > 8.1; See Fig. 2). It should be noted, how- 
ever, that the present models can not show very low /pah 
(< 0.005) for the present galaxies with Ao < 8.1 (i.e., the 
simulated galaxies show low /pah only in their early evolu- 
tion). 

The present study has shown that the radial /pah gradi- 
ents of the present disk galaxies with total masses similar to 
that of the LMG can be only slightly negative (^-^ — 2 x 10~^ 
kpc~^). This almost ffat radial gradients appears to be in- 
consistent with the recent observation by Meixner et al. 
(2010), which has revealed a possible enhancement of /pah 
in the central region of the LMG. The LMG could have ex- 



perienced central starburst events owing to the LMG-SMG- 
Galaxy tidal interaction in the last several Gyr (e.g., Bekki 
& Ghiba 2005), which might enhance the /pah in the central 
regions and consequently enhance the radial /pah gradient. 
The simulated disk galaxies of the present study do not expe- 
rience tidal interaction with other galaxies and thus keep the 
almost ffat /pah gradients. Accordingly, the above appar- 
ent inconsistency could be due largely to differences in the 
evolution of PAH gradients between isolated and interacting 
galaxies (rather than some problems of PAH modeling) . The 
present study thus implies that radial /pah gradients of disk 
galaxies could have some fossil information about the past 
interaction histories with other galaxies. 



4.5 Comparison Avith other models 

The present chemodynamical model incorporated, for the 
ffrst time, the evolution of different chemical components, 
dust formation and destruction, and H2 formation on dust 
grains in a self- consistent manner. Previous simulations, 
however, included H2 formation (but not dust) with different 
models. By assuming that dust abundances are proportional 
to metals (Z) in galaxies, P06 incorporated H2 formation 
on dust grains and H2 destruction by ISRF in SPH hydro- 
dynamical simulations of galaxies for the first time. RK08 
adopted (i) a different model for H2 formation and (ii) a 
more sophisticated ISRF model based on the photoioniza- 
tion code GLOUDY in their SPH simulations of galaxies, 
and discussed global star formation histories of galaxies. 
Gosmological simulations by Gnedin et al. (2009) and Ghris- 
tensen et al. (2012) adopted a model of H2 formation that 
is essentially similar to that by P06. 

Ghemical evolution of variously different components 
(e.g., Mg and Fe) and dust evolution are not included in 
these previous models. The dust-to-gas ratios (D), which 
are key parameters in H2 formation on dust grains, are sim- 
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ply assumed to be proportional to gas-phase metallicities 
(Z) in previous models {Dz — D/Z = constant). As both 
observational and theoretical studies showed that D is not 
simply proportional to Z (i.e., there is large dispersion in D 
for a given Z; Galametz et al. 2011), the adopted assump- 
tion in previous models is neither reasonable nor realistic 
in a strict sense. However, as shown in the present study, 
the time evolution of SFRs and dust properties does not 
depend strongly on whether the evolution of dust is fully 
self-consistently included or not. Thus, the present study 
suggests that adopting a constant Dz in numerical simula- 
tions can be a good approximation for investigating what 
really occurs in real galaxy evolution. 

4.6 Future directions 

Theis & Orlova (2004) incorporated dust-gas interaction in 
hydrodynamical simulations of gas-rich disk galaxies and 
thereby investigated how pressure-free cold dust compo- 
nents, which are assumed to be coupled to the gas by a 
drag force, can influence the gas dynamics of disk galaxies. 
They found that if the dust mass fraction exceeds 2% in 
disk galaxies, then the disk can become destabilized by the 
dust-gas coupling. Although their results suggest dust-gas 
coupling through a drag force could influence the gas dy- 
namics of disk galaxies, such a possibly important effect is 
not included in the present study. Accordingly, it is our in- 
tention to incorporate the dust-gas coupling through a drag 
force in more sophisticated chemodynamical simulations. 

The destruction processes of dust grains by sputtering 
in the reverse shock of a single SN remnant have been in- 
vestigated by hydrodynamical models in detail (e.g., Silvia 
et al. 2010). In comparison with these models on the de- 
tailed physical processes of dust destruction by a single SN, 
the present model for dust destruction could be fairly crude, 
because only the dust destruction time scale is a key param- 
eter for the dust destruction processes. More sophisticated 
parameterization for dust destruction processes will be nec- 
essary so that the time evolution of dust contents can be 
more precisely predicted in our future chemodynamical sim- 
ulations. 



5 CONCLUSIONS 

We have investigated the time evolution of dust, gas, 
and star formation rates in galaxies by using our new 
chemodynamical simulations with a self- consistent model 
for the formation and evolution of dust and molecular 
hydrogen (H2). In this first of a series of papers, we 
have focused particularly on spatial distributions of dust, 
gas-to-dust ratios (D), molecular hydrogen fraction (/h2)5 
gas-phase abundances {Ao =12+log(0/H)), PAH-to-dust 
ratios (/pah), surface mass densities for dust (Sdust), total 
gas (Sgas), and stars (Sstars), and dust-to- stellar mass 
ratios (i^dust = Mdust/Mstar). We have also investigated 
correlations and scaling-relations between these properties 
and the dependences of dust and H2 properties on the halo 
and stellar masses (Mh and Mstar, respectively) of galaxies. 
The principle results are summarized as follows. 

(1) Dust can play a vital role in regulating global star 



formation histories of galaxies, mainly because the time evo- 
lution of H2 formation efficiencies depends strongly on dust 
evolution. Galactic star formation histories thus can depend 
on model parameters for dust formation and evolution, in 
particular, the accretion time scale of dust (race). The ob- 
served D (0.005 ^ 0.01) in luminous disk galaxies can be re- 
produced if Tacc ^ 0.25 Gyr, which is similar to that used in 
previous one-zone chemical evolution models with dust. The 
adopted dust-dependent star formation model is particularly 
important for low-mass disk galaxies with Mh ^ lO^^M©. 

(2) Different local regions of a disk galaxy can have 
significantly different D and Ao owing to their different 
star formation and chemical evolution histories within the 
galaxy. The local regions, however, show a Ao — D corre- 
lation (D oc Ao) and the Ao — D correlation evolves with 
time. The spatial distributions of dust show negative radial 
gradients of D and Ao (i.e., higher D and Ao in inner re- 
gions) in the simulated disk galaxies. 

(3) /h2 rapidly increases owing to more efficient pro- 
duction of dust in the first several Gyr of disk galaxy forma- 
tion. The distributions of /h2 show negative radial gradients 
(i.e., larger /h2 in inner regions) and the gradients evolve 
with time. Local regions with higher D are more likely to 
show higher /h2 . The spatial distributions of H2 gas in disk 
galaxies are more compact than those of dust, and the dis- 
tributions of dust are more compact than those of H i gas. 

(4) Sgas can strongly correlate with Sdust (Sgas oc 
Sdust) in local regions of a disk galaxy. However, Sdust — 
5]h2 correlations are less clearly seen in comparison with 
Sdust — 5]h2 . 5]dust is also well correlated with Spah and 
Sz. i^dust very weakly correlates with Estar in a galaxy such 
that i?dust is higher for lower Sstar. i^dust is higher in early- 
type E/SOs formed by merging in comparison with isolated 
late- type disks. 

(5) If about 5% of dust formed from gaseous ejecta of 
C-rich AGB stars can become PAH dust, then the present 
models can reproduce the observed typical fraction /pah 
[^ 0.03) in luminous disk galaxies. The time evolution of 
/pah is more rapid than that of D^ and /pah does not change 
significantly in the last several Gyr. Weak radial gradients 
of /pah can be seen in the simulated disk galaxies only for 
the first several Gyr evolution, which means that the present 
galaxies are highly likely to have flat or only slightly negative 
/pah gradients. 

(6) Time evolution of SFRs and D is investigated for 
the fiducial model with dust-dependent and metallicity- 
dependent star formation recipes. It is found that there are 
no major differences in the time evolution of SFRs and D 
between the models. These results imply that the evolution 
of SFRs and dust properties are not sensitive to whether a 
dust model is self-consistently included or not for luminous 
galaxies with Mh '^ lO^^M©. It should be stressed, however, 
that the metallicity-dependent star (and H2) formation 
model could possibly underestimate /H25 because it does 
not include the increase of dust-to- metal ratios {Dz) (thus 
the higher probability of conversion from H I to H2) during 
chemical evolution. 

Preliminary results on the dependences of dust and H2 
properties on model parameters (e.g., Mh and Mstar) are 
summarized as follows. 
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(1) The present chemo dynamical model can reproduce 
reasonably well the observed Aq — D relation {D oc Aq), 
though the simulated D for a given Aq is slightly smaller 
than the observed one. Galaxies with larger Mh (Mstar) can 
have higher Aq, D, and /pah owing to more efficient chem- 
ical enrichment (i.e, more efficient dust production) in more 
massive galaxies. 

(2) The final (i.e., present) values of i^dust depend 
strongly on Mh (and Mstar) such that galaxies with larger 
Mh (and larger Mstar) have lower i^dust- Early- type galax- 
ies transformed from late-type disks via galaxy merging can 
have lower i?dust in comparison with late- type disks. Thus, 
the present study suggests that morphological transforma- 
tion of galaxies can cause significant evolution of i?dust • 

(3) Disk galaxies with larger Mh are more likely to have 
higher /h2 • The final /h2 can be very low in low-mass disk 
galaxies with Mh for Mh < lO^^M©, because star formation 
is more severely suppressed in these low-mass disks (i.e., a 
smaller amount of dust is produced and used for H2 forma- 
tion). 

(4) Disk galaxies with larger Mdust have larger Mgas 
(= Mhi + Mhs) and the simulated Mdust — Mgas relation is 
consistent reasonably well with the observed one (Mgas oc 
^dust)- The Mdust — Mh2 relation derived for the simulated 
disk galaxies is slightly steeper than the observed relation 
(Mh^ocM^VsI). 
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Table Al. A summary of parameters of subgrid physics for 
additional test models. 



model 


fchl 


A^ISRFF _ 


Otuiayi 


fiducial 


1.0 


1.0 


1.4 


Tl 


0.5 


1.0 


1.4 


T2 


2.0 


1.0 


1.4 


T3 


1.0 


1.0 


0.7 


T4 


1.0 


1.0 


2.1 


T5 


1.0 


1.5 


1.4 


T6 


1.0 


2.0 


1.4 



^ The parameter A^h is used for determining column densities 
of neutral and molecular hydrogen: Ni i = / !" * ni idr ^ 
2fch^i,i^i for neutral hydrogen. 

The parameter /cisfr is used for determining the local vol- 
ume VisRF (oc ^fgj^pe|) of the surrounding region of a gas 
particle. 
^ The maximum time step width in units of 10^ yr. 
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APPENDIX A: NUMERICAL TESTS FOR 
SUBGRID PHYSICS 

In the present study, H2 formation on dust grains and H2 
destruction by ISRF are estimated for each SPH particle at 
a given time step by using the smoothing length (hi) and the 
gravitational softening length (cg). For example, we estimate 
the column density of H i gas (which is necessary to estimate 
H2 formation) as follows: 



/khhi 



ni,idr ^ 2khni,ihi, 



(Al) 



where kh is set to be 1 in all models described in the main 
text. Accordingly we need to demonstrate that the present 
results do not depend on the resolution of simulations (i.e., 

/Ch). 

ISRF for a gas particle is estimated for the local volume 
VisRF that is defined as; 



V^iSRF = A:isRFeg, 



(A2) 



where eg is the gravitational softening length for the gas 
particle and /cisrf is set to be 1 for all models described in 
the main text. Such local estimation of ISRF implies that 
stellar radiation from stars that are located outside eg (for 
A:isFR = 1) is completely ignored. Although, the radiation 
density would not change significantly even if a larger vol- 
ume is used, we need to demonstrate that the present results 
do not depend on /cisrf- Hopkins et al. (2011) have discussed 
a similar resolution problem related to the estimation of the 
total amount of radiation from massive stars in their simu- 
lations. Also, the present results could possibly depend on 
the maximum time step width ((5tmax), though timescales of 
key physical processes (e.g., Tdust) are much longer than the 
adopted (5tmax for models described in the main text. 

We thus investigate the models in which model param- 
eters are the same as those adopted in the fiducial model 
except /ch, /cisRF, and ^tmax- The parameter values for these 
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Figure Al. The same as Fig. 11 but for the seven comparative models (including the fiducial model): comparison between different k-^ 
(left), 5fmax (middle), and /cisrf (right). 



test models (Tl— T6 as well as the fiducial model) are sum- 
marized in Table Al. Key results on the time evolution of 
D and /hs, which are shown in Fig. Al, are summarized 
as follows. First, the time evolution and final values of D 
and /h2 do not depend on k\, for 0.25 ^ k\, ^ 2.0. This 
result means that the present results are insensitive to the 
method to estimate column densities of neutral and molecu- 
lar hydrogen. Second, the final fu^ value in the model with 
(^tmax = 0.7 Myr (a finer time step width) is 34% larger 
than that in the fiducial model, though there are no signifi- 
cant differences in D between the two. The final /h2 in the 
model with ^tmax = 2.1 Myr is 25% smaller than that in the 
fiducial model. 



Gas densities can become higher in the model with 
smaller ^tmax so that the probability of conversion from neu- 
tral to molecular hydrogen can become higher. As a result 
of this, /h2 can become larger. The derived dependences of 
/h2 on 5tmax imply that fn^ can be better predicted by us- 
ing models with smaller Stma^. However, a larger amount of 
calculation time is required for completing numerical simu- 
lations with smaller ^tmax- Thus, if one has a limited amount 
of allocated time for performing numerical simulations (on 
GPU clusters or supercomputers), one would need to com- 
promise on the accuracy of the predicted /h2 • 



Fig. Al also shows that the final /h2 and D do not de- 
pend so strongly on /cisrf: /h2 is slightly higher in the mod- 
els with larger /cisrf, probably because ISRF (i.e., radiation 
density) can be somewhat weaker for a larger volume (i.e., 
more efficient /h2 formation due to weaker ISRF). These 
results imply that the present results are not so sensitive to 
the adopted method to estimate ISRF (i.e., non-inclusion of 
the contribution of distant stellar particles to ISRF). 



APPENDIX B: 
DISKS 



SF HISTORIES OF LOW-MASS 



Fig. Bl shows the time evolution of SFRs in the low-mass 
model with Mh = lO^^M© for the H- and H2-dependent 
SF recipes. Although the fiducial model shows only slight 
differences in SFRs between the two different SF recipes, 
this model shows more remarkable differences in the early 
SF histories. The model with the H— dependent SF shows 
systematically higher SFR for T = 1 ^-^ 2 Gyr in compar- 
ison with that with the H2— dependent SF. This suppres- 
sion of SF in early phase of the disk formation in the H2- 
dependent SF model is due largely to the lower H2 fraction 
in the low-mass disk model. Such suppression can be seen 
in other low-mass models with the H2— dependent SF recipe 
(Mh ^ lO^^M©). These results imply that the adoption of 
H2— dependent SF recipe is important for SF histories of 
low-mass disk galaxies. Severe suppression of SF has been 
already reported by K12 in which H2-regulated SF models 
are adopted. 



APPENDIX C: LOWER Fa^ IN GALAXIES 
WITH LOWER TOTAL MASSES 

In the present model, H2 formation efficiencies depend 
strongly on dust abundances (D). Therefore, it is expected 
that more massive galaxies, where chemical enrichment pro- 
ceeds more efficiently, can have larger D and thus larger 
/h2- Fig- CI shows the time evolution of D and /h2 in the 
low-mass disk model with Mh = 10^° M© and the high-mass 
one with Mh = lO^^M©. Clearly, both final D and /h2 are 
higher in the high-mass model, and the evolution of /h2 is 
more rapid in the high-mass model. In the present H2 forma- 
tion model, the final /h2 in disks with Mh ^ 10^° M© can be 
very low (< 0.05), which has important implications of gas 
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Figure Bl. The time evolution of SFRs in the two models with 
H-dependent (blue) and H2 dependent (red) SF recipes. These 
results are for the low-mass disk models with M^ = IQ-'^'-'M© and 
A = 0.038. 
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Figure CI. The same as Fig. 11 but for the low-mass halo 
model with M^ = 10-*^^ M© (blue) and the high-mass model with 
Mh = IQi^M© (red). 



contents in dwarf galaxies. We will discuss such implications 
in detail in our forthcoming papers. 



